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ABSTRACT 

The  aerodynamic  performance  of  turbomachines  is  limited  by  flutter,  choking,  and  stall,  the  last  of  which 
being  a  topic  of  this  report.  Viscous  flow  considerations  are  not  only  essential  in  determining  loss  coefficients, 
but  also  in  ascertaining  loadings  and  in  providing  a  fundamental  understanding  of  the  flow.  However,  because 
of  the  complexity'  of  flows  in  turbomachines  flow  problems  of  model  character  only  are  accessible  to  uumerical 
treatment.  The  focus  of  this  report  is  on  steady,  incompressible,  viscous.  2-d  flows  through  a  single  infinite 
row  of  equidistant  blades.  In  order  to  compute  these  cascade  flows  Prandtl's  boundary  layer  concept  enhanced 
by  a  strong  interaction  model  is  applied  to  the  Navier-Stokes  equations.  Since  this  approach  enables  the 
integration  of  boundary  layer  equations  in  the  presence  of  backflow.  it  is  an  attractive  alternative  to  much 
more  expensive  Navier-Stokes  solvers. 

Cebeci  and  collaborators,  who  developed  such  a  viscous-inviscid  interaction  method  for  single  airfoils. 
modified  and  extended  their  code  to  the  cascade  problem.  The  present  contribution  consists  of  a  detailed 
description  of  the  theory*  and  of  conducting  a  few  sample  calculations  of  test-character.  Flows  past  staggered 
cascades,  using  the  ten  percent  thick  airfoils  of  the  NACA-65-series  compressor  blade  family  as  cross  sections, 
were  studied  over  a  range  of  Reynolds  numbers  and  air  inlet  angles.  Difficulties  were  encountered  in  cases  of 
low  Reynolds  number  and  highly  loaded  blades.  The  apparent  success  of  the  single  airfoil  code  and  its  close 
agreement  with  experiments  could  not  yet  be  extended  to  the  cascade  problem,  partly  because  of  the  more 
demanding  airfoil  shapes  used  in  turbomachines  and  partly  because  of  the  inadequacy  of  the  turbulence 
model  at  low  Reynolds  numbers.  Improvements  concerning  both  numerics  aud  physics  are  necessary  to  turn 
the  current  cascade  code  into  a  robust  and  efficient  tool  to  analyse  cascade  flows  of  practical  interest. 
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1.  INTRODUCTION  AND  GENERAL  CONCEPTS 

StalL  the  phenomenon  at  which  lift  decreases  with  an  increase  in  angle  of  incidence,  is  one  of  the  limiting 
factors  in  the  efficient  operation  of  turbomachines.  The  flow  regime  at  stall  is  characterized  by  flow  separation 
leading  to  regions  of  reversed  flow,  which  do  not  only  determine  the  loss  behaviour  but  cause  significant 
changes  in  loading  and  overall  lift  force.  Flow  separation  is  also  present  in  another  flow  regime.  Le.  bubble 
transition,  which  is  a  low  Reynolds  number  flow  phenomenon,  comprising  laminar  separation,  transition,  and 
turbulent  reattachment.  Because  of  economical  requirements  on  turbomachines,  which  result  in  higher  blade 
loadings,  the  aerodynamic  community  is  forced  to  provide  a  better  understanding  of  viscous  effects  in  both 
attached  and  separated  flow.  The  increasing  costs  of  wind  tunnel  experiments  and  the  decreasing  computing 
costs  encourage  the  development  of  new  numerical  tools.  Solutions  of  the  full  Navier-Stokes  equations  are 
still  reserved  for  supercomputers,  and  on  the  other  hand,  classical  boundary  layer  codes  are  restricted  to 
attached  flows.  Therefore,  engineers  have  developed  procedures,  which  do  the  job  of  a  Navier-Stokes  solver  at 
a  lower  expense  of  computing  time  and  storage.  Such  procedures  are  based  on  the  viscous-inviscid  interaction 
concept  which  will  be  used  in  the  present  analysis. 

What  are  the  ingredients  of  a  successful  viscous-inviscid  interaction  method?  Following  Prandtl's  bound- 
ary layer  concept  the  computational  domain  is  divided  into  an  in  viscid  region  and  a  thin  layer,  close  to  the 
blade's  surface,  where  viscous  effects  cannot  be  neglected.  The  computational  scheme  solves  separately  for 
the  flows  in  the  viscous  and  inviscid  domains,  which  was  already  a  feature  of  classical  boundary  layer  codes. 
but  the  essential  point  is  the  mechanism  how  inviscid  and  viscous  flows  are  matched.  While  classical  bound- 
ary layer  methods  suffer  from  a  pronounced  hierarchy  of  viscous  and  inviscid  regions,  interaction  models 
are  supposed  to  exhibit  the  ability  of  mutual  and  immediate  transfer  of  information  between  viscous  and 
inviscid  regions.  In  summary,  interaction  methods  must  include  an  inviscid  flow  solution,  a  boundary  layer 
method,  able  of  handling  both  attached  and  separated  flow,  and  some  interaction  model. 

At  this  point  we  have  to  concede  that  our  numerical  possibilities  are  not  even  close  to  a  solution  of  the 
flow  inside  a  turbomachine.  Unsteadiness,  compressibility,  highly  3-d  character,  interaction  of  shocks  and 
wakes,  and  complex  flow  boundaries,  as  found  in  turbomachines.  move  a  solution  of  the  real  problem  far 
beyond  today's  potential.  Because  of  the  need  to  come  up  with  some  results,  model  problems  are  introduced. 
which  retain  only  a  limited  number  of  the  original  features.  Regarding  geometry  and  decomposition  of  the 
3-d  problem  we  adopt  Wu's  ideas  to  split  the  actual  flow  into  two  separate,  but  interacting  2-d  flows,  namely 
a  blade-to-blade  flow  and  a  meridional-through  flow.  The  blade-to-blade  or  cascade  plane  is  defined  as  a 
surface  of  revolution  formed  by  rotating  a  streamline  of  the  meridional  plane  about  the  machine's  axis. 
Such  a  projection  of  a  single  blade  row.  and  our  considerations  are  confined  to  single  blade  rows,  reveals  an 
infinite  row  of  equidistant,  similar,  airfoil-shaped  bodies,  called  cascade.  Since  we  do  not  attempt  a  solution 
of  the  3-d  problem,  the  cascade  problem  will  be  isolated  from  the  meridional-through  flow.  Further  we 
do  not  allow  any  upstream  disturbances,  as  generated  by  the  wakes  of  upstream  blade  rows.  The  current 
calculations  will  also  neglect  the  effect  of  viscous  wakes  emanating  from  the  trailing  edges  of  the  considered 
blades.  There  is  no  doubt  that  the  flow  in  a  turbomachine  is  highly  unsteady  and  even  the  artificial  cascade 
problem  might  experience  unsteadiness  due  to  a  distorted  upstream  flow  field.  However,  because  our  model 
presumes  steady  flow,  our  ability  to  predict  stalled  flows  will  be  somewhat  restricted.  Our  limitation  to 
incompressible  flows  is  mainly  related  to  the  availability  of  efficient  inviscid  cascade  codes.  This  assumption 
dispenses  us  from  taking  shocks,  shock  reflections,  and  shock  interactions  into  account.  Finally  the  boundary 
layer  approach  imposes  the  restriction  that  viscous  effects  remain  confined  to  a  thin  layer.  Obviously,  this 
assumption  is  violated  in  case  of  extensively  separated  flows.  In  view  of  all  these  drawbacks  what  advantages 
are  being  offered  by  inviscid-viscous  interaction  methods  and  why  is  the  interaction  method  superior  to 
classical  boundary  layer  methods  and  purely  inviscid  flow  solvers?  The  following  levels  of  information  can  be 
associated  with  these  procedures:  Inviscid  flow  solvers  yield  a  pressure  distribution  and  the  overall  lift  force, 
which  can  be  expected  to  be  accurate  only  in  case  of  attached  and  high  Reynolds  number  flows.  Boundary 
layer  methods  give  additional  information  on  skin  friction  distribution  and  overall  drag  force,  but  they  fail  if 
the  flow  separates.  Interaction  methods  provide  pressure  and  skin  friction  distributions,  lift  and  drag  forces 
for  attached  and  mildly  separated  flows. 

The  following  sections  will  shed  some  light  on  the  individual  parts  of  this  approach.  The  inviscid  flow 
equations  are  being  solved  by  a  standard  panel  method,  which,  mathematically  speaking,  provides  a  solution 
of  Laplace's  equation  subject  to  Neumann  type  boundary  conditions.  Singularity  methods,  such  as  the  panel 


method,  attempt  a  solution  by  an  appropriate  superposition  of  elementary  solutions,  named  singularities,  so 
that  the  boundary  condition  is  satisfied.  In  two-dimensional  flows  only  two  types  of  singularities  are  required, 
namely  sources  ( and  sinks)  to  model  dbplacement  flows,  and  vortices  to  model  lift.  There  remains  the  choice 
where  and  how  to  distribute  these  singularities.  Smith  and  Hess  [13]  proposed  to  place  them  on  the  boundary 
surfaces  which  led  to  the  first  and  still  very  successful  panel  codes.  The  first  order  approximation,  which  we 
are  going  to  use.  employs  a  representation  of  the  airfoil  by  straight-line-elements,  along  which  singularity 
strengths  are  assumed  constant.  The  airfoil  code  of  Smith  and  Hess  was  extended  to  cascades  by  Giesing  [9j, 
whose  approach  is  virtually  the  same  as  the  current  one,  except  that  Giesing's  method  is  derived  and  coded 
in  complex  variables,  while  the  current  one  is  done  in  real  variables.  Inviscid  methods  intended  to  interact 
with  viscous  flow  solvers  require  some  modifications.  The  commonly  used  flow  tangency  condition  has  to  be 
replaced  by  the  prescription  of  a  wall  transpiration  velocity  to  make  allowance  for  the  displacement  effect 
due  to  the  existence  of  boundary  layers. 

The  second  building  block  of  the  scheme  is  the  viscous  flow  solution,  being  achieved  by  a  finite  dif- 
ference method.  Finite  difference  methods  are  numerical  tools  to  solve  partial  differential  equations  by 
discretizing  the  unknown  functions  at  a  limited  number  of  gridpoints  and  subsequent  rewriting  of  the  differ- 
ential equations  in  form  of  algebraic  equations,  each  of  which  applying  to  an  elementary  subdomain.  The 
algorithm  depends  on  the  propagation  of  pertubations,  such  that  elliptic  equations  require  expansive  field 
solutions,  while  parabolic  equations  permit  marching  procedures.  Hence  engineers  try  to  avoid  the  elliptic 
Navier-Stokes  equations  and  to  focus  on  the  parabolic  boundary  layer  equations.  An  order  of  magnitude 
estimation,  making  use  of  the  relation  that  the  boundary  layer  thickness  is  much  smaller  than  the  streamwise 
length  coordinate,  justifies  the  reduction  of  the  elliptic  Navier-Stokes  equations  to  the  parabolic  boundary 
layer  equations.  However,  no  gain  without  penalties:  Only  flow  fields  which  obey  the  boundary  layer  ap- 
proximation ought  to  be  treated,  computations  will  break  down  at  the  point  of  zero  skin  friction  due  to  the 
Goldstein  singularity,  and  the  downstream  marching  procedure  will  face  numerical  instabilities  in  regions  of 
reversed  flow.  The  first  limitation  is  inherent  to  the  approach.  However,  the  Goldstein  singularity  is  related 
to  the  boundary  conditions.  Catherall  and  Mangier  [3]  discovered  that  the  boundary  layer  equations  can  be 
integrated  through  the  point  of  separation,  if  displacement  thickness  or  wall  shear  stress  is  prescribed  in- 
stead of  the  external  velocity.  The  third  problem  is  tricky  because  of  both  up-  and  downstream  propagating 
pertubations.  which  revert  the  problem  locally  to  an  elliptic  one.  The  most  common  remedy,  the  so-cailed 
FLARE  approximation  [19],  drops  the  streamwise  convective  term  in  regions  of  backflow.  thus  getting  rid 
of  the  numerical  instabilities,  but  violating  the  physics  whenever  backflow  velocities  exceed  small  values. 
After  having  discussed  the  basic  advantages  and  disadvantages  of  boundary  layer  methods,  we  continue  with 
information  on  the  currently  employed  method.  Our  viscous  flow  solution  features  Keller's  box  method  [14j. 
which  is  implicit  and  second  order  accurate.  In  regions  where  separation  is  unlikely  to  occur  the  boundary 
layer  equations  are  recast  in  terms  of  Falkner-Skau  variables  and  solved  subject  to  prescribed  external  veloci- 
ties. This  direct  method,  which  is  applied  in  the  immediate  neighbourhood  of  the  leading  edge  only,  exhibits 
the  classical  hierarchical  structure  between  viscous  and  inviscid  flow  regions,  indicating  that  the  outer  flow 
imposes  the  pressure  on  the  viscous  layer.  Soon  after  the  leading  edge  or  at  the  pressure  peak,  at  the  latest, 
the  scheme  switches  from  the  direct  to  the  interactive  mode,  which  ensures  the  mutual  interaction  of  inner 
and  outer  flow.  While  direct  methods  specify  the  pressure  and  inverse  methods  the  displacement  thickness, 
neither  of  these  is  prescribed  in  an  interactive  boundary  condition.  Both  external  velocity  and  dbplacement 
thickness  are  treated  as  unknowns  and  the  boundary  condition  specifies  a  relation  between  the  unknown 
streamfunction  at  the  boundary  layer  edge,  external  velocity  and  displacement  thickness.  The  concept  of 
nearly  constant  boundary  layer  thickness,  as  provided  by  Falkner-Skan  variables,  has  to  be  abandoned.  The 
boundary  layer  equations  are  now  written  in  terms  of  physical  transformed  coordinates  and  the  Mechul  func- 
tion method  by  Cebeci  accounts  for  the  status  of  the  external  velocity  as  an  unknown  variable.  Besides  the 
crucial  boundary  couditiou  at  the  outer  edge  of  the  viscous  layer  both  methods  employ  the  no-slip  condition 
at  the  wall  and  the  direct  method  requires  initial  conditions,  which  are  generated  by  similarity  solutions. 

In  order  to  compute  turbulent  flows  knowledge  about  Reynolds  stresses,  specifically  turbulent  shear 
stresses,  is  necessary.  Unfortunately  our  understanding  of  turbulence  is  still  quite  inadequate.  Therefore. 
we  must  use  simple  models  based  on  hypotheses.  The  currently  used  model,  the  eddy  viscosity  concept  by 
Cebeci  and  Smith  [4j.  relates  the  turbulent  shear  stresses  to  mean  flow  properties  on  an  empirical  basis.  This 
algebraic  zero-equation  closure  allows  to  retain  the  same  form  of  the  differential  equations  and  to  apply  the 


same  numerical  procedures  to  both  laminar  and  turbulent  flows.  In  order  to  know  where  to  start  turbulent 
flow  computations  prediction  of  transition  is  necessary.  Since  again  no  exact  methods  are  available,  the  best 
choice  is  to  take  experimental  results.  The  attempt  to  predict  transition  is  still  in  its  infancy,  resulting  in 
rather  limited  applicability.  Michel's  method,  which  is  used  here,  is  based  upon  an  empirical  data  correlation, 
narrowing  considerably  the  range  of  applicability,  while  linear  stability  theories  with  empirical  assumptions, 
like  the  e9-method  offer  some  universality.  Probably,  transition  and  turbulence  modelling  represent  the 
weakest  members  in  the  computational  scheme. 

The  interaction  term  describes  the  way  how  the  inner  viscous  and  outer  inviscid  regions  communicate 
with  each  other.  The  transmitter  of  any  information  is  obviously  the  edge  boundary  condition,  which  in  case 
of  a  weak  interaction  works  out  as  follows:  The  result  of  the  inviscid  flow  computation,  which  is  the  pressure 
distribution,  is  imposed  as  boundary  condition  on  the  viscous  flow  computation.  In  return  the  result  of  the 
viscous  flow  solution,  which  is  the  displacement  thickness,  determines  the  effective  shape  of  the  body  which 
is  input  to  the  inviscid  flow  solver.  Basically,  there  are  two  possibilities  to  account  for  the  displacement  effect 
of  a  boundary  layer  on  the  outer  inviscid  flow:  Either  the  blades  are  thickened  by  adding  the  displacement 
thickness  to  the  blade  coordinates  or  the  flow  tangency  condition  on  the  blade  is  substituted  by  a  condition 
which  prescribes  a  blowing  velocity.  Both  ways  can  be  incorporated  either  in  a  weak  or  strong  interaction.  A 
weak  interaction  will  be  sufficient  if  the  viscous  effect  on  the  pressure  is  smalL  permitting  the  above  described, 
hierarchical  manner  of  solving  the  complete  flow  problem.  However,  there  are  regions,  like  the  trailing  edge 
and  separated  flow,  where  viscous  displacement  might  cause  more  severe  changes  in  the  local  pressure  field 
so  that  alternating  treatment  of  external  velocity  and  displacement  thickness  as  in-  and  output  appears  to 
be  too  rigid  a  mechanism  for  the  interaction  of  viscous  and  inviscid  flow  regimes.  These  sensitive  areas  call 
for  equivalent  and  simultaneous  treatment  of  displacement  thickness  and  external  velocity.  Schemes,  which 
meet  these  standards  and  hence  do  not  possess  a  definite  hierarchy,  are  named  strong  interaction  models. 
Only  they  are  capable  of  calculating  boundary  layers  with  flow  separation.  The  approach  taken  here  follows 
closely  the  ideas  of  Veldman  (22),  who  initially  introduced  the  concept  of  strong  interaction. 

The  global  organization  of  the  scheme  (see  figure  1  on  page  4)  shows  an  iterative  process  in  which 
alternately  inviscid  and  viscous  flow  equations  are  being  solved.  Results  are  considered  as  converged  if  the 
external  velocity  distributions  of  viscous  and  inviscid  flow  solutions  do  coincide  and  if  the  results  of  the 
viscous  flow  solutions  do  remain  steady  over  consecutive  iteration  cycles.  Relaxation  is  a  numerical  tool 
either  to  ensure  or  to  accelerate  convergence.  Viscous  flow  results  are  updated  by  relaxation  just  before  they 
are  transferred  to  the  inviscid  flow  solver  to  serve  there  as  boundary  condition.  The  actually  used  relaxation 
formula  was  proposed  by  Kwon  and  Pletcher  [15]  and  allows  overtaxation. 

The  present  report  expands  upon  the  work  of  Cebeci  et  aL  [8|,  who  developed  a  viscous-inviscid 
interaction  method  for  incompressible  cascade  flows.  The  objective  of  the  present  work  is  to  provide  a 
complete  description  of  the  theory  and  to  conduct  a  few  sample  calculations  in  order  to  begin  the  evaluation 
procedure.  In  chapter  2  a  description  of  the  panel  method  is  given  which  is  used  to  compute  incompressible 
inviscid  flows  through  arbitrarily  staggered  cascades  of  arbitrarily  shaped  blade  sections.  The  solution  of 
the  boundary  layer  equations  subject  to  an  imposed  pressure  as  well  as  to  an  interactive  edge  condition  is 
discussed  in  chapter  3.  which,  additionally,  reviews  the  used  interaction  and  turbulence  models.  Flow  results 
in  form  of  pressure,  skin  friction,  displacement  thickness,  and  velocity  profiles  are  presented  and  analysed 
in  chapter  4  for  the  NACA-65-series  compressor  blade  family.  The  cascade  configuration,  namely  a  stagger 
angle  of  45  degrees  and  a  blade  solidity  of  (3.78.  was  selected  in  consideration  of  available  experimental  data. 
Calculations  were  performed  for  Reynolds  numbers  of  245  thousand  and  6  million  and  for  air  inlet  angles 
ranging  from  40  to  62  degrees.  The  report  closes  with  a  brief  summary  and  recommendations  for  future 
work.  The  reader  is  alerted  to  the  fact  that  the  present  report  is  of  an  interim  nature  covering  the  period 
September  1985  through  August  1986. 
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Fig.  1.  Global  organization  of  viscous-inviscid  interaction  method 


2.  INVISCID  FLOW  METHOD 

Since  the  outer  flow  is  assumed  to  be  incompressible,  a  solution  of  Laplace's  equation  will  provide  in  viscid 
flow  results.  Laplace's  equation  will  be  subject  to  a  Neumann  type  boundary  condition,  which  prescribes 
a  blowing  velocity  on  the  surface  of  the  blades,  resulting  from  the  viscous  displacement  effect.  In  addition 
Kutta's  condition,  postulating  zero  load  at  the  trailing  edge,  is  imposed  to  establish  a  unique  solution. 
The  numerical  implementation  requires  a  couple  of  approximations  which  are  associated  with  the  term 
panel  method.  The  following  sections  will  cover  a  description  of  singularity  methods  in  potential  flows,  the 
definition  of  the  inviscid  cascade  problem,  the  concepts  of  influence  coefficients  and  basic  flows,  and  the 
derivation  of  farfield  cascade  relations  enabling  the  statement  of  Kutta's  condition  for  given  inlet  angle  of 
attack  and  the  calculation  of  inlet  and  exit  flow  conditions. 

2.1  Singularity  methods  In  potential  flows 

Irrotational  flow  fields  ( curl  V  =  o)  permit  that  the  velocity  vector  may  be  expressed  as  the  gradient  of  a 
scalar  function,  called  the  velocity  potential  *, 

V  =  gradQ. 

This  definition  together  with  the  principle  of  conservation  of  mass  for  incompressible  flows  [die  V  =  0)  yields 
Laplace's  equation, 

div[grad$)  =  0        or        A$  ■  0 

which  governs  the  steady  flow  field  of  incompressible  inviscid  fluids.  The  commonly  used  flow  tangency 
condition  is  replaced  here  by  a  boundary  condition  that  specifies  a  wall  transpiration  velocity 

-—  =  IV.W*)        on  the  surface  of  blades 
on 

Laplace's  equation  is  a  linear  homogeneous  second  order  partial  differential  equation,  which  is  classified  as 
elliptic,  demanding  a  field  solution.  Two  concepts  neve  proven  successful  in  the  exact  solution  of  Laplace's 
equation,  the  method  of  separation  of  variables  and  singularity  methods.  The  latter  ones  achieve  a  solution 
by  an  appropriate  superposition  of  elementary  solutions,  called  singularities.  This  technique,  underlying  all 
panel  codes,  is  based  on  Green's  theorem  and  the  principle  of  superposition,  which  applies  to  any  linear 
partial  differential  equation 

If  A*i  =  0.  A$.2  =  0 A<&„  =  0 A$:v  =  0 

then  A*  =  0  with  $  =  (7t#i  +  (72*8  + ...  +  Gn$n  +  .. .  +  C.v*  v 

where  C\,  C2, ....  (7„, . . . ,  and  Cy  denote  arbitrary  constants.  Green's  theorem,  which  is  specific  to  Laplace's 
equation,  states  that  the  flow  past  arbitrarily  shaped  bodies  may  be  represented  by  surface  singularity 
distributions.  This  theorem  simplifies  the  solution  procedure  considerably,  because  the  flow  at  a  field  point 
can  be  computed  by  summing  up  the  effects  of  all  singularities  instead  of  all  field  points.  For  2-d  flows  the 
singularities  to  be  used  are  sources  and  vortices,  because  the  effect  of  displacement  and  thickness,  respectively. 
can  be  handled  by  distributions  of  sources  and  sinks,  while  the  effect  of  lift  is  modelled  by  a  distribution  of 
vortices.  The  name  singularity  is  attached  to  sources  and  vortices,  because  they  do  satisfy  Laplace's  equation 
everywhere  except  at  a  singular  point. 

Since  the  overall  flow  will  be  built  up  by  simple  flows,  we  review  first  the  velocity  potentials  of  basic 
flows  which  constitute  the  elements  of  our  flow  synthesis.  The  velocity  potential  of  uniform  parallel  flow 
that  is  inclined  with  respect  to  the  positive  x-axis  by  an  angle  a  is  given  by 

■p[  x.  y )  =  17«o  ( x  cos  o  +  y  sin  a ) 

where  U^  denotes  the  magnitude  of  the  constant  velocity  vector.  Radial  streamlines  aud  concentric  circular 
equipotentials  characterize  the  flow  field  of  a  single  source 
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where  <r  denotes  the  source  strength  and  the  source  is  located  at  (0,0).  Interchanging  streamlines  and 
equipotentials  of  a  single  source  produces  the  flow  field  of  a  single  vortex 

<p{*<  9)  ■  -  =£  arctan  j 

where  a  positive  7  denotes  the  strength  of  a  clockwise  rotating  vortex,  located  at  (0,0).  Our  principal  interest 
is  in  cascade  flows  so  that  our  further  considerations  will  deal  with  infinite  rows  of  singularities.  The  velocity 
potential  of  an  infinite  row  of  equal  and  equidistant  sources  at  the  points  (0, 0),  (0,  ±r),  (0,  ±2r),  (0,  ±3r), . . . 
can  be  expressed  in  terms  of  hyperbolic  and  trigonometric  functions 

"%    fc=+cc  *"        L  *  '  " 

A  similar  process  leads  to  the  velocity  potential  of  an  infinite  row  of  equal  and  equidistant  vortices  at  the 
points  (0,0),  (0,±r),  (0,±2r),  (0,±3r),... 

*=-Ko                    .                            /  tan  —  \ 
<fi{*,y)  =  - 57     E     arctan  V— —  ■  - ^-arctanl  £7-) 

r 

2.2  Cascade  flows 

Recalling  that  a  cascade  is  defined  as  an  infinite  row  of  equidistant  similar  airfoil  shaped  bodies,  we  can 
conclude  that  the  flow  field  past  a  cascade,  which  is  inclined  to  a  uniform  steady  onset  flow,  consists  of 
identical  subflowflelds.  each  of  which  is  enclosed  by  two  streamlines  being  displaced  at  the  cascade's  spacing. 
Identical  subflowflelds  obviously  can  be  caused  only  by  identical  singularity  distributions  on  each  body,  so 
that  the  above  defined  infinite  series  will  replace  the  expressions  for  single  source  and  vortex  as  used  in  the. 
airfoil  formulae.  On  the  basis  of  Laplace's  equation's  linearity  we  decompose  the  total  velocity  potential  $ 
into  that  of  uniform  onset  flow  and  the  so-called  disturbance  potential  55,  which  accounts  for  the  disturbance 
by  the  cascade  of  an  originally  uniform  flow  field 

^(j*.  y)  =  Vyz  (•»*  cos  a  +  y  sin  a)  +  p(ar,  y) 

where  U.^.  and  a  are  magnitude  and  direction  (with  respect  to  positive  x-axis)  of  the  vector  mean  velocity. 
(Note:  The  vector  mean  velocity  is  defined  as  the  vector  average  of  inlet  and  exit  velocity  vector).  The 
disturbance  velocity  potential  of  a  cascade  .  whose  blades  are  arranged  along  the  y-axis.  can  be  written  as 


So    * 


/  ta.n-(y- y0(s))  \  } 

-  —arctan    '— )  ds 

2»  Vtanh-U-JoUn/J 


where  the  integration  is  carried  out  along  the  surface  of  the  center  blade,  (jo«//o)  denotes  the  coordinates 
of  the  center  blade,  and  r  is  the  vertical  spacing  between  adjacent  blades.  Since  this  integration  cauuot  be 
performed  analytically  for  arbitrarily  shaped  blades,  the  coutiuuos  singularity  distributions  are  discretized 
in  the  following  fashion: 

1.  Blades  wUl  be  represented  by  straight  line  panels. 

2.  The  boundary  condition  will  be  satisfied  at  the  midpoints  of  panels. 

3.  Source  strengths  are  assumed  constant  along  each  panel,  but  vary  from  panel  to  panel.  Vortex  strengths 
are  constant  over  the  blade. 

This  unequal  status  of  source  and  vortex  strengths  is  due  to  the  requirement  that  the  number  of  unknowns 
must  balance  the  number  of  equations.  If  the  blade  is  approximated  by  N  panels.  N  boundary  conditions 
applied  at  the  midpoints  of  the  X  panels,  plus  Kutta's  condition  at  the  trailing  edge  produce  equations 
for  iV  +  1  unknowns.    Hence  the  N  unknown  source  strenghts  are  determined  such  that  the  boundary 
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conditions  are  satisfied,  while  there  is  only  one  equation  left  to  adjust  one  vortex  strength.  Introducing  these 
approximations  to  the  disturbance  potential  we  get 

tp(x,y)=    C  t£    /    In    -f  cosh — (x-Xj-tcoaBj)  -  cos — (y-y,-*sin0,)J    rf* 
Milt 
"    +///a  /  tan  -( y-  y,  -  «  sin  fj ) 


-~E      /    arctan    £ d* 

>=1    ?m  VtMh -(#-#/-•  COily)/ 


where    <r,  denotes  the  source  strength  of  the  j-th  panel 

7  the  vortex  strength,  which  is  constant  for  all  panels 

[xj,  yj)  are  the  coordinates  of  the  midpoint  of  the  j-th  panel 

lj  and  8 }  denote  length  of  the  j-th  panel  and  the  angle  ( positive  counterclockwise)  between  surface 

direction  of  the  j-th  panel  (positive  clockwise)  and  positive  x-axis. 

2.3  Concept  of  influence  coefficients 

Eventually  we  are  concerned  with  velocities  rather  than  the  potential  which  is  merely  a  vehicle  to  give 
information  on  the  velocity  field.  Remembering  that  the  velocity  vector  is  the  gradient  of  the  potential  we 
can  apply  this  operator  in  any  convenient  reference  coordinate  system  including  local  coordinate  systems 
attached  to  panels.  Now  there  are  two  velocity  components  we  need  to  know:  Firstly,  normal  velocities 
at  panel-midpoints  are  required  to  satisfy  the  boundary  conditions  and  secondly,  tangential  velocites  are 
necessary  to  compute  the  pressure  distribution  on  the  blades.  Both  velocities  will  be  computed  by  separately 
taking  into  account  the  influence  of  each  panel  resulting  in  2.V  contributions,  one  half  induced  by  source 
distributions,  the  other  half  induced  by  vortex  distributions.  This  procedure,  called  concept  of  influence 
coefficients,  constitutes  the  numerical  core  of  the  panel  method.  Influence  coefficients  are  hence  defined 
as  velocities  induced  by  a  unit  singularity  distribution  of  one  paneL  Let  «,-  and  tf,  be  outer  normal  and 
clockwise  tangential  directions  at  the  midpoint  (/,. yj  of  the  i-th  panel  then  influence  coefficients  due  to 
source  distributions  take  the  form: 
Aij  represents  the  outer  normal  velocity  at  (/,.  y.)  induced  by  a  unit  source  distribution  on  the  j-th  panel 

1      f      d    f     f  1  /         2*  2r  \  ] 

Aij  ~  T~    j      ]T~  \  lu    ^\  cosu  ~H *i  ~ *j ~  *  cos 9 j )  -  cos  — ( u> -  !Jj-  *  sin  0, )  J    >  ds 

'-,j/2       'I  > 

=  -  sin(  $i  -  9j)  VZiJ  +  cos(  Bt  -  Bj)  VyiJ         for  i  £  j 

Aij  =  -  for  i  -  j 

Bij  represents  the  clockwise  tangential  velocity  at  (*•,#«)  due  to  a  unit  source  distribution  on  the  j-th  panel 

+<j/2        fro  o  1) 

Bij  *=  ~    /     —  <  In    -(  cosh  y{  x, ■-  xj  -  s  cos  Bj )  -  cos  y{  U\  -  !J}  -  *  sin  B} )  J    >  ds 

■■     '      -0/2         '    I  ' 

=  cos(  Bi  -  Bj)  YXij  -  sin|  Bt  -  Bj)  Vyij        for  i  £  j 
Bij  -  0  for  i  =  ; 

with 

2  vjxj-Xj)          2  zjyt-Yj) 
i  cosh — * cos — — 

VXil  =  —In 


Vg.j  —  —  I  arctanAi  +  arctanA?  -  arctanAjj 


Vjfj  and  Vy.j  denote  x-  and  y-components  ( in  a  local  coordinate  system,  whose  x-  and  y-axes  point  in  the 
clockwise  tangential  and  outer  normal  directions  of  the  j-th  panel)  of  the  velocity  at  (*,-,  yt ).  induced  by  a 
unit  source  distribution  of  the  j-th  paneL  Coordinates  of  boundary  points  (A",,  Yj)  and  midpoints  (/>,#,) 
are  related  by 

*i  -  \i  *i  +  A'y+i )  -  Xj  +  -L  cos  9j 

Potential  flow  solutions  frequently  encounter  troubles  because  of  the  multivalued  arctangent.  The  following 
two  provisions  make  sure  that  the  arctangent,  and  subsequently  the  induced  velocity,  do  not  jump  in  value, 
where  it  is  not  supposed  to  do  so:  All  arctangents  must  be  evaluated  in  the  range  —  tc  to  +x  and  V9iJ  is 
calculated  by  adding  up  the  effect  of  a  single  airfoil,  given  by  the  first  term,  and  the  effect  of  a  cascade 
without  the  center  blade,  represented  by  the  second  and  third  term.  A  minor  difficulty,  the  computation 
of  velocities  at  the  midpoint  of  the  inducing  panel,  can  be  overcome  by  Cauckys  principal  value  technique, 
which  permits  integration  of  certain  singular  integrands. 

The  influence  coefficients  due  to  a  vortex  distribution  can  be  related  to  those  due  to  a  source  distribution. 
Based  on  the  fact  that  velocity  fields  of  vortex  and  source  are  different  by  a  clockwise  rotation  of  90  degrees 
only,  the  normal  velocity  due  to  a  vortex  equals  the  negative  tangential  velocity  due  to  a  source,  and  the 
tangential  velocity  due  to  a  vortex  equals  the  normal  velocity  due  to  a  source 


iVORTEX  _       n5' 


SOURCE 


nVORTEX  _  . SOURCE 

Bij  ~       Atj 

(Unless  indicated  otherwise  influence  coefficients  will  be  those  due  to  a  source  distribution.) 

2.4  Concept  of  basic  flows 

Equipped  with  the  concept  of  influence  coefficients  we  can  now  deal  with  the  satisfaction  of  the  boundary 
condition,  which  is  to  determine  source  and  vortex  strengths  such  that  velocities  of  given  magnitude  and 
direction  be  induced  at  panel  midpoints: 

=  ( v  „,,),  for  «  =  1 .\ 


an, 
-V  .v 

JO) 


J2  AeJ<r      -7  E  Bij ■  —  ( ]'u-t)t  +  r/odsintfj  sina  -  cosfl, ;coso) 


for  /  =  1 .V 


To  match  the  ,V  +  1  unknowns  ( .V  source  and  1  vortex  strength)  by  an  equal  number  of  equations  the  system  is 
complemented  by  the  empirical  Kutta's  condition.  The  resulting  system  could  be  solved  in  a  straightforward 
fashion  by  a  linear  equation  algorithm:  however,  traditionally  the  flow  is  split  into  three  basic  flows,  thus 
providing  a  very  efficient  solution  for  different  angles  of  attack  in  case  of  no  wall  transpiration.  Basic  flows 
naturally  divide  into  displacement  and  lift  generating  flows  as  follows: 
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1.  Zero  angle  of  attack  flow  (or  s  0  deg,  7  =  0,  wall  transpiration): 


jv 

E 


/si        y 


(0)  _  (^n'f)t 


cos  a 


+  Una  sin  9i        for  t  ■  1, . . . ,  N 


2.  Ninety  degrees  angle  of  attack  flow  (a  m  90  deg,  7  ■  0,  no  wall  transpiration): 


s 


/=i 


for  1  ■  1, . . . ,  JV 


3.  Purely  circulatory  flow  {U^  =  0,  7  =  1.  no  wall  transpiration): 


E  4*4 

/si 


(7)  =   E  *,/ 
/si 


for  1  ■  1 iV 


Each  of  the  above  equations  stands  for  a  system  of  N  equations  with  unknown  source  strengths  only.  Since 
all  systems  show  up  with  the  same  coefficient  matrix,  the  three  sets  of  source  strengths  can  be  obtained  by 
one  execution  of  a  linear  equation  solver  with  three  different  right  hand  sides.  Although  the  concept  of  basic 
flows  is  primarily  a  numerical  tool  the  following  physical  interpretation  can  be  given:  The  first  basic  flow 
represents  a  pure  displacement  flow,  whereby  blades  are  inclined  at  zero  angle  of  attack  ( angle  with  respect 
to  global  x-axis.  not  blade  chordline).  The  boundary  condition  of  the  first  basic  flow  consists  in  a  prescriptiou 
of  a  wall  transpiration  velocity.  The  second  basic  flow  is  also  a  pure  displacement  flow,  but  at  an  angle  of 
attack  of  ninety  degrees  and  subject  to  a  no-penetration  condition.  The  third  basic  flow  represents  a  purely 
circulatory  flow.  Since  a  constant  vortex  distribution  alone  would  not  obey  the  flow  tangency  conditiou. 
normal  velocities  induced  by  this  constant  vortex  distribution  are  balanced  by  an  extra  source  distribution. 
In  summary:  Each  of  the  basic  flows  satisfies  some  kind  of  boundary  condition  and  each  of  them  violates 
Kutta's  condition.  An  appropriate  linear  combination  of  all  three  basic  flows  will  yield  flow  at  the  required 
angle  of  attack  subject  to  Kutta's  coudition: 

Total  flow        =         "Zero  angle  of  attack  flow"  *  cos  a 

+  "Ninety  degrees  angle  of  attack  flows"  *  sin  a 

+  "Purely  circulator^"  flow**  *  7 

The  vortex  strength  7  is  still  unknown  in  the  above  equation,  but  provided  it  is  known,  total  normal  and 
tangential  velocities  at  the  midpoint  of  the  i-th  panel  take  the  following  form 


Vn{Xi,yi)  ■ 


1 1(  *■•&-)  = 


diij 


V  V 

•  =  f/co  sin(  cr  -  9i)  +   E  Aiy<rj  -  7  E  BtJ 

7=1  ;=l 

.V  .V 

=  r/co  cos(a  -  9i)+   E  Bijffj  +  7  E  Atj 

/Sl  J  =  L 


with  <Tj  =  <r(0)  cos  a  +  <r(9° 


sin  a  +  <r.    7. 


The  deviation  of  the  outer  normal  velocities  Vn  from  the  prescribed 
wall  transpiration  velocities  Vwt  is  a  measure  of  the  accuracy  of  the  linear  equation  solver. 

2.5  Farfleld  cascade  relations 

For  farfleld  considerations  the  flow  field  past  a  cascade  can  be  simulated  by  an  infinite  row  of  discrete  vortices. 
Such  an  approach  can  be  justified  because  the  displacement-disturbances  decay  with  increasing  up-  and 
downstream  distance  and  because  there  remains  a  finite  effect  due  to  circulation  in  the  farfleld.  Using  the 
formerly  defined  velocity  potential  of  an  infinite  row  of  equal  and  equidistant  vortices,  this  simplified  model 
produces  information  about  the  horizontal  and  vertical  disturbance  velocities  in  the  farfleld  of  a  cascade 


Hi  =  it  2  =     Urn 

SS—S.30 


2 


tan 


-  — -arctau 


.   si 
tauh  — 


tan 


—  —  arctanl   irjr 

2'T  V  tanh  — 


where  the  small  letters  u  and  v  denote  the  components  of  disturbance  velocities,  and  the  subscripts  1  and 
2  indicate  streamwise  positions  at  up-  and  downstream  infinity.  Magnitudes  and  angles  of  total  inlet  and 
exit  velocities  can  be  obtained  by  superposition  of  the  above  disturbance  flowfield  and  a  uniform  paralell 
fiowfield,  which  is  inclined  by  an  angle  of  a 


Abs(  Inlet  velocity):        \Jv\  +  Vj  -  y/(  UZ  cos  o  +  Ux )2  +  ( UZ  sin  a  +  Vx  )2 


-\M£) 


+  — 


...,,,  ,  V,       t^siii  or +  I7(2r) 

Air  mlet  angle:  tane»i  =  — -  =  — « 

"i  "co  cos  Of 


.   Abs( Outlet  velocity):       >Jv*  +  V?  =  Ju^  +  (£)*"  u*>™aT 

V2       U^sma-T/{2r) 

Air  outlet  angle:  tan  a  a  =  — -  ■  — '■ — — ' 

Ui  Urx,  cos  a 

where  the  capital  letters  U  and  V  designate  the  components  of  total  velocities  and  Uoe  is  the  magnitude  of 
the  vector  mean  velocity.  (For  an  illustration  of  angles  and  velocity  components  see  figure  2  on  page  12). 

The  calculation  of  forces  requires  either  a  pressure  integration  along  the  blade  or  the  application  of  the 
momentum  equation  in  integral  form  to  an  appropriate  control  surface.  The  latter  method  is  to  be  preferred 
here  because  it  involves  fewer  numerical  steps  and  is  said  to  be  more  accurate  for  small  numbers  of  panels. 
Assuming  steady  2-d  flows  of  incompressible  fluids  the  momentum  equation  can  be  written  according  to 


f  m  (F„  Fy)T  =  -pi  V[Y1\)  ds-  I  pn  d$ 


where  C  denotes  the  control  surface  consisting  of  two  vertical  planes  in  the  up-  and  downstream  farfield  and 
two  similar  streamlines,  which  are  arranged  such  that  they  enclose  a  single  blade.  Making  use  of  Bernoulli's 
equation  magnitude  and  direction  of  the  resultant  gross  force  per  unit  span  become 


F=^F]  +  Fi=pU^r        with    U^  =  <Ju?  +  ((\\  +  V2)/2)2 
e  =  arctan(Fy/Fx)  =  a  +  90  degrees 

The  gross  force  turns  out  to  act  at  a  right  angle  relative  to  vector  mean  velocity,  while  no  force  is  exerted 
in  the  direction  of  vector  mean  velocity.  Lift  is  defined  as  the  force  perpendicular  to  vector  mean  velocity  . 
Normalizing  this  force  by  a  reference  area  and  the  dynamic  pressure,  based  either  on  vector  mean  or  inlet 
velocity,  gives  the  dimeusionless  lift  coefficient 

n  F  2U*X 


Cref(pU%f/2)         CrefU%f 

2.6  Kutta's  condition 

Kutta's  condition  assures  that  flows  on  the  upper  and  lower  surfaces  merge  smoothly  at  the  trailing  edge  by 
postulating  equal  pressures  and  subsequent  equal  velocities  there.  Since  the  panel  method  does  not  permit 
the  evaluation  of  velocities  at  boundary  points.  Kutta's  condition  is  satisfied  in  an  approximate  fashion  at 
the  coutrol  points  of  the  first  and  last  panel 

V**,  =  -VfV 

This  condition  determines  the  vortex  strength  7.  whose  knowledge  is  necessary  to  synthesize  the  combined 
flow.  Hence  the  scheme  makes  use  of  Kutta's  condition  after  completing  the  basic  flow  solution  and  before 
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starting  the  combined  flow  computation.  Cascade  flows  offer  several  input  possibilities,  two  of  which  are 
being  considered  here.  For  a  given  average  onset  flow  angle  Kutta's  condition  is  a  linear  equation  in  the 
unknown  vortex  strength  7 

ffoo[cos(a-0i)  +  cos(a-0iv)]+   E  [(Biy  +  BNi){cosa<r\0)  +  sin  a  <r\90)  )^ 

r    s  S  i 

=  -7     E  (B1>+B^7)+   E  (-Aw  +  4v>) 


For  given  inlet  air  angle  two  unknowns,  average  onset  flow  angle  a  and  vortex  strength  7,  remain  to  be 
determined,  resulting  in  the  need  for  another  relation.  Kutta's  condition  and  the  equation  for  the  inlet  air 
angle  constitute  a  linear  system  in  the  two  unknowns  -7/  cos  a  and  tan  a,  which  can  be  solved  by  any  suitable 
procedure 


y  ,  y  y 


*  >=1  J  L  jssi  *  >=1 

=  -^oo[cos^+cos^)]-    E  (Bij  +  B/v.,V,0) 


7 

cos  a 


S/  •      7 

C^oo  tan  a  +  ■—*■  =  ^  tan  oti 

2r    cos  o 

One  possibility  to  increase  accuracy,  especially  if  there  are  lower  limits  on  the  panel  length  at  the  trailing 
edge,  is  extrapolation  of  velocities  at  the  trailing  edge.  This  approach  avoids  the  problem  regarding  singular 
behaviour  of  velocities  at  panel  boundary  points,  and  allows  to  impose  Kutta's  condition  exactly  at  the 
trailing  edge.  In  case  of  linear  extrapolation  Kutta's  condition  can  be  written  as 

«i  +  h  iy  +  »;v-i 

where  the  indices  1  and  X  refer  to  the  first  and  last  paueL  with  numbering  starting  at  the  trailing  edge  ou 
the  lower  surface  and  ending  there  ou  the  upper  surface. 

If  the  displacement  effect  due  to  the  boundary  layer  is  accounted  for  ( nonzero  wall  transpiration  velocity), 
it  is  advisable  to  evaluate  the  pressure  distribution  on  the  displacement  surface  itself.  This  step  has  to 
be  taken  because  if  the  wall  transpiration  velocity  exceeds  small  values,  the  assumption  of  110  pressure 
gradient  accross  the  ficticious  inviscid  flow  region  between  blade  and  displacement  surface  ceases  to  be 
valid.  The  displacement  surface  is  therefore  the  location,  where  total  velocities  are  to  be  computed  and 
where  Kutta's  condition  Is  to  be  met.  While  the  scheme  up  to  and  including  the  calculation  of  basic  flows 
remains  unchanged,  the  scheme  then  proceeds  with  the  computation  of  two  more  sets  of  influence  coefficients. 
These  so-called  off-body  influence  coefficients  represent  normal  and  tangential  velocities  on  the  displacement 
surface  induced  by  source  and  vortex  distibutions  of  unit  strength,  residing  on  the  original  blade  surface. 
The  scheme  makes  use  of  these  coefficients  in  determining  circulation  by  imposing  Kutta's  condition  off- 
body  at  the  trailing  edge  and  computing  total  velocities.  On  the  assumption  of  a  rigorous  (inviscid)  theory 
the  displacement  surface  will  not  be  penetrated  by  fluid  particles.  However,  because  the  computation  of 
the  blowing  velocities  is  based  upon  a  flat  plate  model,  total  velocities  on  the  displacement  surface  have  a 
tangential  velocity  component  and  a  nonzero  normal  velocity  component,  which  is  smaller  by  an  order  of 
magnitude.  The  magnitude  of  the  total  velocities  multiplied  by  the  sign  of  the  tangential  component  are 
transferred  to  the  viscous  flow  solver. 


n 


Upstream 
Farfield 


Downstream 
Farfield 


Fig.  2.  Cascade  geometry  and  its  vector  diagrams 


LEGEND 


j 

j-l/2 
H 


l 

K/2 

ki/2 

f 

V 

i 

i 

^ J 

hj/2 

/ 

!                   ,  hjyfe 

S                  '                  XT 

V 

1 ► 

O     known 
D     unknown 

A     center   for  momentum 
equation 

V     center  for   equations 
containing 
q-derivatives   only 


i-l  i-l/2  i 

Fig.  3.  Net  rectangle  for  finite  difference  approximations 
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Fig.  4.  Visualization  of  the  blowing  velocity  concept 
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8.  VISCOUS  FLOW  METHOD 

In  1904  L.  Prandtl  gave  birth  to  boundary  layer  theory  which  closed  the  gap  between  the  exact  science  of 
hydrodynamics  and  the  more  engineering-orientated  hydraulics.  Prandtl  showed  that  the  flow  field  past  a 
body  can  be  divided  into  two  regions  : 

1.  a  thin  layer  close  to  the  surface  of  the  body  and  the  ensuing  wake,  where  viscous  effects  have  to  be 
accounted  for,  and 

2.  the  remaining  region,  where  inertia  terms  are  dominant  over  viscous  terms  so  that  the  latter  may  be 
neglected. 

These  simplifications  turn  the  elliptic  Navier-Stokes  equations  into  the  parabolic  boundary  layer  equations, 
which  restrain  the  upstream  propagation  of  disturbances.  Numerically,  this  change  of  characteristics  trans- 
forms the  algorithm  from  a  field  procedure  to  a  marching  method,  which  integrates  the  boundary  layer 
equations  for  given  initial  conditions  step  by  step  proceeding  in  the  downstream  direction.  The  differential 
equations  are  complemented  by  boundary  conditions,  whose  choice  significantly  influences  the  applicability 
of  each  scheme: 

1.  Direct  boundary  layer  methods  employ  the  so-called  no-slip  condition,  requiring  zero  normal  and  zero 
tangential  velocity  at  the  surface,  and  a  prescription  of  the  external  velocity  at  the  edge  of  the  boundary 
layer.  Such  procedures  encounter  the  Goldstein  singularity  [10]  at  the  point  of  separation,  leading  to  a 
numerical  breakdown. 

2.  Boundary  layer  equations  can  be  integrated  through  the  point  of  separation,  if  the  above  edge  boundary 
condition  is  replaced  by  a  prescription  of  either  displacement  thickness  or  wall  shear  stress,  while  the 
no-slip  condition  is  kept  unchanged.  However,  such  an  inverse  boundary  layer  method  coupled  with  an 
inverse  inviscid  flow  solution  suffers  from  very  poor  convergence. 

3.  Interactive  boundary  layer  methods  utilize  an  edge  boundary  condition  which  possesses  the  elliptic 
character  of  the  outer  flow.  This  so-called  interaction  law,  which  prescribes  a  linear  combination  of  dis- 
placement thickness  and  external  velocity,  and  the  no-slip  condition  provide  the  required  three  boundary 
conditions  for  this  scheme. 

The  numerical  approach  features  a  finite  difference  method,  which  recasts  the  continuity  and  momentum 
equations  as  a  system  of  linear  algebraic  equations.  Kellers  box  method  [14]  approximates  derivatives  by 
central  differences  and  linearizes  the  nonlinear  terms  around  the  solution  of  the  previous  iteration.  To 
stabilize  the  numerical  algorithm  in  the  presence  of  backflow  the  FLARE-approximation  [19]  is  adopted  in 
which  the  streamwise  convection  term  is  dropped  in  regions  of  backflow.  The  resulting  set  of  difference 
equations  is  solved  by  Keller's  block  elimination  method  which  takes  advantage  of  the  sparsely  occupied 
matrix 

In  the  following  descriptions  of  both  direct  and  interactive  boundary  layer  methods,  a  derivation  of 
the  interaction  law.  and  minimum  information  on  the  turbulence  model  will  be  given.  The  section  "Direct 
Boundary  Layer  Method"  will  elucidate  the  more  basic  ingredients,  like  discretization  of  the  flow  field 
and  Newton's  linearization,  while  the  section  on  interactive  boundary  layer  methods  will  emphasize  its 
pecularities. 

3.1  Direct  boundary  layer  method 

Direct  methods  are  applicable  in  the  absence  of  flow  reversals,  which  limits  their  use  in  the  current  approach 
to  the  immediate  neighbourhood  of  the  leading  edge.  They  are  supposed  to  generate  initial  conditions  at 
the  stagnation  point  and  to  allow  an  efficient  integration  of  the  boundary  layer  equations  around  the  leading 
edge.  Generally  speaking,  direct  methods  appeal  to  the  weak  interaction  problem  in  which  viscous  effects  on 
the  pressure  are  required  to  remain  small.  Although  the  paper  addresses  the  strong  interaction  problem,  all 
primary  features  of  the  finite  difference  method  may  as  well  be  demonstrated  for  the  direct  method.  To  begin 
with  we  consider  steady  2-d  flows  of  incompressible  fluids,  described  in  a  curvilinear  coordinate  system  with 
/  directing  along  the  blade's  surface  and  y  denoting  the  outer  normal  direction.  The  velocity  components  u 
and  v  shall  be  determined  such  that  they  satisfy  anywhere  in  the  flowfield  the  equations 

dn       dv 

dx       dy 

dn         dn  diie  d  . ,  dn . 

u-r-  +  i>  —  -  ue-j-  +  t/—[b—) 
dx         dy  ax         dy     dy 
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and  at  the  boundaries  of  the  flowfield  the  conditions 

y=0:        «(x,0)=0,         v(x,0)  =  0 
y=ye:      «(*♦»«)*««(*) 

where  b  denotes  l  +  vt/v.  The  equations  are  referred  to  as  boundary  layer  or  thin  shear  layer  equations,  the 
conditions  as  no-slip  and  edge  condition.  In  spite  of  the  box  method's  requirement  to  reduce  the  governing 
equations  to  a  system  of  first-order  equations,  we  introduce  for  the  present  a  streamfunction  (u  =  dii'/dy 
and  v  =  -di'/dx),  which  naturally  decreases  the  number  of  dependent  variables.  Since  streamfunctions  by 
definition  satisfy  the  continuity  equation,  only  the  momentum  equation  is  left,  which  becomes  a  third  order 
equation  in  d' 

d^    d2il>       d&  d34>  _      due         d  {JP*x 

dy  dxdy      dx   dy2         e  dx         dy     dy2 

Application  of  the  Falkner-Skan  transformation  contributes  to  the  removal  of  the  singularity  at  /  =  0  and 
maintains  a  nearly  constant  boundary  layer  thickness  in  terms  of  the  new  variable.  In  contrast  to  similarity 
concepts  the  dependency  on  two  independent  variables,  x  and  q.  is  kept.  The  Falkner-Skan  transformation 
scales  the  normal  coordinate  y  and  the  streamfunction  if  with  reference  to  local  external  velocity  and  local 
streamwise  coordinate 


"v£» 


With  prime  denoting  differentiation  with  respect  to  i\  the  momentum  equation  and  boundary  conditions 
become 

(&/T+^/r+m[i-(/ri  =  /(/'|^-r^) 

//  =  0:        /'(*,0)    =0.      /(z,0)  =  0 

where  m  —  (x/ue)(due/dx)  is  a  dimensionless  pressure  gradient  parameter. 

The  box  method  involves  four  steps:  reduction  to  a  first  order  system,  discretization  in  terms  of  central 
differences  and  two-point  averages,  linearization,  and  solution  of  the  resulting  linear  system.  The  first  step  is 
accomplished  by  the  introduction  of  two  additional  dependent  variables.  U  and  V,  which  converts  the  third 
order  momentum  equation  into  a  system  of  three  first  order  equations 

f'  =  U 
U'  -  V 

The  system  is  complemented  by  the  boundary  conditions,  which  read  in  terms  of  the  new  variables 

t)  sQ:       U[x.O)   =0.     /(x,0)  =  0 
i}  =  He:      U[x.i\e)  =  l 

The  second  step  deals  with  the  discretization  of  the  continuous  flow  problem  by  considering  instead  of  the 
continuos  functions  /.  U .  and  V  a  set  of  discrete  values  of  these  flow  properties.  Let  the  solution  domain 
0  <  x  <  xi,     0  <  //  <  i) j  be  covered  by  a  rectangular  mesh 

xy  =  0.  j,  =  /,_!  +  kj      with     2<  i  <  I 

Ii  =  0.  i]j  =  ilj-i  +  hj     with     2  <  j  <  J     and     r\j  =  t\e 
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Any  flow  property,  whether  it  is  a  function  or  its  derivative,  will  be  expressed  in  terms  of  nodal  values  and 
coordinates  of  the  network.  Triplets  of  flow  properties,  consisting  of  streamfunction  and  its  first  and  second 
derivative  with  respect  to  q,  are  assigned  to  each  node  of  the  network 

{/(*,,*,),  u[*i,m),  v(xhnj))  =  {/;,  tf«,  v}} 

Because  of  the  parabolic  behaviour  of  the  boundary  layer  equations  the  solution  at  a  certain  streamwise 
position,  say  /,,  depends  solely  on  the  solutions  of  upstream  positions,  say  ?i_i,  x ,_a, : •. .,  while  no  down- 
stream influence  has  to  be  considered.  The  overall  solution  works  in  a  step-by-step  fashion,  propagating  in 
the  downstream  direction.  The  use  of  first  order  equations  and  the  technique  of  central  differences  allow  to 
reduce  the  domain  of  dependence  from  all  upstream  x-stations  to  the  immediate  preceding  one.  Hence  one 
step  of  the  solution  procedure  writes  the  governing  equations  for  a  column  of  net  rectangles  (called  boxes) 
residing  in  the  subdomain 

**— i  <  «p  <  *t     and    0  <  q  <  qj 

and  solves  subsequently  for  the  nodal  values  at  the  downstream  face  of  the  rectangular  shaped  subdomain. 
The  x-station  being  currently  solved  holds  therefore  the  superscript  T,  while  the  superscript  "/- 1"  denotes 
the  known  flow  properties  of  the  adjacent  upstream  location.  As  indicated  by  the  term  central  differences 
the  equations  are  satisfied  midway  between  the  nodes: 

The  two  ordinary  equations  are  centered  about  (x,,  n}-\ji ).  and 
the  partial  equation  is  centered  about  [xt-i/tilf-if*)* 
The  problem  can  now  be  rewritten  in  terms  of  finite  differences  (see  figure  3  on  page  12) 
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Boundary  conditions  undergo  the  same  discretization  process  and  take  the  form 


U{  =  0, 


ft-Q 


One  might  question  whether  this  represents  a  solvable  system  in  which  the  number  of  unknowns  matches 
the  number  of  equations.  The  subdomain  under  consideration  consists  of  J  — 1  net  rectangles,  the  flow 
quantities  of  each  being  related  by  a  momentum  equation.  Two  auxiliary  relations,  which  link  dependent 
variables  to  their  q -derivatives,  refer  to  the  downstream  face  of  each  net  rectaugle.  Recalling  the  three 
boundary  conditions  we  have  a  total  of  ZJ  equations,  which  face  J  triplets  of  unknowns,  each  of  the  triplets 
corresponding  to  anode  of  the  downstream  side  of  the  rectangular  subdomain.  Unfortunately,  the  unknowns 
appear  in  nonlinear  combinations,  so  that  the  solution  involves  an  iterative  procedure.  The  variables  are 
linearized  around  their  values  of  the  preceding  iteration  (and  in  case  of  the  first  iteration  around  the  solution 
of  the  adjacent  upstream  station) 


J  j         ~  Jj 

u'-K  =  u)-K-1 

T,'«.K    T -I.K— 1 
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where  k  denotes  the  iteration  counter.  If  these  Newton  iterates  are  introduced  and  subsequently  terms  that 
are  quadratic  in  [6fjK,  6Uj,K,  6Vj'K)  are  dropped,  a  linear  system  in  the  unknowns  6fj *,  6UjK,  6Vj'K  is 
obtained 

tup  -  nrgi  -  ^(*v;-*  +  sv;*)  =  u^r1  -  *}**  +  *y »$5J 

(*i)J"«f  +  MjVlft  +  (*)/'*<*  +  (•4)</'**/jf1  +  («s#*  «7^  +  («t)J"«£i  =  («•*)'/* 
where  the  coefficients  abbreviate  the  following  expressions 
ft*:*-1 


"r        ft.    ~  '     >~V2   'J- 1/2    "*"    *  j-l/2Jj-l/2    "   Jj-l/2  S-l/2  ' 


^       Vi/a;;-i/«  +  "m 

Together  with  the  boundary  couditious  this  system  is  repeatedly  solved  until  the  Newton  iterates  6f''K,  6U'.'K, 

6Vj'K  are  small  enough  to  be  neglected.  As  frequently  experienced  in  discretized  boundary  value  problems 
the  resulting  set  of  linearized  finite  difference  equations  displays  a  block-tridiagonal  structure  which  permits 
a  fast  algorithm  being  applied  to  the  linear  equation  solution.  Block-tridiagonal  matrices  are  composed  of 
submatrices.  named  blocks,  of  which  only  those  residing  on  main  and  both  adjacent  diagonals  have  nonzero 
entries.  Equations  and  boundary  conditions  of  the  discretized  boundary  layer  problem  require  a  special 
arrangement  to  exhibit  this  block-tridiagonal  structure.  The  system  consists  of  J  sets,  each  of  which  include 
three  equations: 

a.  First  set  of  equations  (index  j  =  1) 

1.  Wall  boundary  coudition  prescribing  no  penetration 

2.  Wall  boundary  condition  prescribing  zero  tangential  velocity 

3.  Second  auxiliary  relation  (linking  U'  and  V")  of  the  first  box 

b.  Intermediate  sets  of  equations  (index  2  <  j  <  J  —  l) 

1.  First  auxiliary  relation  (linking  /'  and  U)  of  the  (J-l)-st  box 

2.  Momentum  equation  of  the  (/- l)-st  box 

3.  Second  auxiliary  relation  (linking  V  aud  V)  of  the  J-th  box 
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c.  Last  set  of  equations  (index  j  m  J) 

1.  First  auxiliary  relation  (linking  /'  and  17)  of  the  last  box 

2.  Momentum  equation  of  the  last  box 

3.  Edge  boundary  condition. 

Boxes  are  numbered,  the  first  referring  to  that  at  the  wall,  the  last  referring  to  that  at  the  boundary 
layer  edge.  With  [Sf)  =  (f^.iOf*,^}'  and  {r«:*}  =  {(rl)iJK,{r1t)ijKAri)i-K}T  denoting  vectors  of 
the  unknown  Newton  iterates  and  the  known  right  hand  sides,  respectively,  the  linearized  system  can  be 
expressed  according  to 


i<t*i 
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where  the  blocks  are  3  x  3-matrices  being  defined  as 
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and  where  the  right  hand  sides  are  obtained  from 


10  0 

0      1  0 

[0    -1     -ft2/2j 


UV*1  = 


*8ir 


-A//2 
1 


0 

>i)y" 
0 


fl'.K-l 


•  <.K-l 


•  I.K-l 


nij   -  fj-i    -  f j      +"j<-'j-i/2 

as  given  by 


;ra)^  = 


momentum  equation 


'j 

ii.K 


=  u 


«.A^l 


'4./C—  I 


7  ■        UJ-l        +llJ+^j 


«.K-l 

+  1/2 


(n)i=0,  (r2)i-*  =  0.         (r3)^  =  0 


for  2  <  ;  <  J 
for  2  <  j  <  J 
for  1  <  j  <  J  -  1 


The  variety  of  iudices  explains  itself  by  the  2-dimensionality  and  uouliuear  character  of  the  problem: 

1.  The  superscript  T  indicates  the  location  along  the  blade  surface. 

2.  The  subscript  "/*  denotes  the  location  perpendicular  to  the  blade  surface. 

3.  The  superscript  "k"  designates  the  iteration  counter. 

The  unknowns  are  the  Newton  iterates  of  streamfunction  and  its  first  and  second  derivative  at  the  /'-th 
streamwise  position  of  the  K'-th  iteration  all  across  the  boundary  layer.  The  above  system  of  linear  equations 
can  be  solved  most  efficiently  by  Keller's  block  elimination  method,  which  consists  primarily  of  two  types  of 
recursions.  A  forward  step  eliminates  the  lower  diagonal  of  submatrices.  wheras  the  backward  step  solves 
the  remaining  system  from  bottom  to  top. 
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Note  on  discretization:  The  averaging  of  all  terms  is  carried  through  such  that  the  number  of  generated 
terms  is  a  minimum,  which  turns  out  to  be  trivial  for  single  quantities,  but  has  some  impact  on  the  product 
terms.  Furthermore,  averaging  makes  no  difference  whether  acting  on  dependent  variables  or  dependent 
parameters  like  "m"  and  "&".  Let  a  and  d  be  any  dependent  flow  quantities,  then  discretization  is  performed 
on  the  basis  of  the  following  rules: 

1.  Evaluation  of  single  terms:  a;-i/2  m  z(a)  +  a)-i  +  <1y~1  +  a/-i) 

/da\'-V*      a' -a*'1 

2.  Evaluation  of  differential  quotients:     (  —J  = 

/da\  _  fly-  Oy-! 

V^iy/y-i/a  hj 

3.  Evaluation  of  product  terms:  (a<0y-i/2  =  o-[(<M9/+  (a<0y-i]  =  *[*&  +  aj-idj-i) 

rather  than:  (<*<0y-i/2  —  ji aj^j  +  ajdj-\  +  fly-irfy  +  a/.idy.i) 

Although  the  second  way  of  evaluating  the  product  term  is  just  as  valid  as  the  first  one,  the  more  straightfor- 
ward fashion  of  the  first  formula  is  preferred  throughout  the  method.  Care  is  advised  in  the  reverse  process, 
so     <* j- 1/2  rfy-1/2  7*  ( <*<0j- 1/2     because  of  the  above  agreement. 

3.2  Interactive  boundary  layer  method 

The  integration  of  boundary  layer  equations  subject  to  an  imposed  pressure  distribution  terminates  at  the 
point  of  zero  skin  friction.  Moreover,  direct  methods  are  restricted  to  weakly  interacting  regions,  where 
viscous  disturbances  to  the  inviscid  flow  remain  smalL  Strong  interaction  arises  for  example  from  boundary 
layer  separation  or  the  rapid  flow  acceleration  downstream  of  the  trailing  edge,  both  of  which  cause  sub- 
stantial changes  in  the  external  velocity  distribution  due  to  viscous  displacement.  Interactive  methods  are 
designed  to  permit  integration  through  the  point  of  zero  skin  friction  and  to  account  for  strong  interaction. 
The  boundary  layer  approach  is  maintained  in  strongly  interacting  regions,  but  accommodations  apply  to 
the  status  of  external  velocity  and  edge  boundary  condition.  In  contrast  to  direct  and  inverse  methods  both 
external  velocity  and  displacement  thickness  are  treated  here  as  unknown  quantities,  reflecting  the  elliptic 
character  of  the  outer  flow.  This  introduces  apparently  one  additional  unknown  into  the  viscous  flow  prob- 
lem, whose  solution  can  be  achieved  either  by  an  eigenvalue  method  or  the  Mechul  function  method,  the 
latter  of  which  being  preferred  here.  The  edge  boundary  condition  of  the  direct  problem  is  supplemented 
by  the  so-called  interactive  boundary  condition,  which  relates  the  unknown  external  velocity  with  its  invis- 
cid and  "displacement-perturbation-related"  contributions.  Boundary  layer  equations  in  the  following  form 
constitute  a  system  in  the  unknown  functions  u[x,  /y),  i>(x,y),  and  ue(x,  y) 

dn       dv 

dx       ay 

dn  dn  dne  d  . , dn . 

n  - —  +  i>  —  =  ne (-  v — (0  — ) 

ox         dy  dx  dy     dy 

0=  -r^ 

The  system  consists  of  continuity  equation,  x-momentum  equation,  and  the  seemingly  unnecessary  y- 
momeutum  equation,  wherein  pressure  in  the  momentum  equations  has  been  expressed  in  terms  of  the 
external  velocity.  The  Mechul  function  approach  assumes  that  the  external  velocity  be  a  function  of  two 
arguments,  resulting  in  the  need  for  the  trivial  y-momeutum  equation.  The  reason  for  considering  »e(jr.  y) 
rather  than  \%t\  x)  is  of  a  purely  numerical  nature,  i.e..  such  a  provision  allows  an  easy  setup  of  finite  difference 
equations,  avoiding  nonlinear  eigenvalue  techniques.  The  governing  equatious  are  complemented  by  proper 
boundary  conditions:  The  velocity  components  u  and  t>  are  required  to  satisfy  the  no-slip  couditiou  at  the 
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surface  of  the  blade  and  the  horizontal  component  must  merge  smoothly  into  the  outer  flow  at  the  boundary 
layer  edge 

y=0:        «(x,0)  =0,     t'(*,0)  =  0 


y=lfc 


d* 


u(x,ye)  m  ue(x,ye),     ue(x,y€)  »  uet(x)  +  -  /-jt(M*) — 

t  J  at  x  - 
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with  ««/(/)  denoting  the  inviscid  velocity  distribution  and  the  second  term,  named  Hilbert  integral,  approx- 
imating the  perturbation  velocity  due  to  viscous  effects. 

Interactive  methods  can  be  employed  in  both  attached  and  separated  flow  regions  .while  direct  methods 
cannot  do  so  because  of  their  breakdown  related  to  the  Goldstein  singularity  nor  inverse  methods  because  of 
their  poor  convergence  rates.  Since  interactive  methods  offer  such  a  favourable  behaviour,  they  are  used  in  the 
current  scheme  for  the  main  parts  of  the  blade.  Only  the  stagnation  point  singularity  prohibits  an  exclusive 
use  of  the  interactive  method.  The  steps  which  turn  the  partial  differential  equations  of  the  interactive 
problem  into  a  linear  system  of  algebraic  equations  resemble  those  of  the  direct  method  so  that  only  major 
steps  and  their  results  will  be  repeated  here.  After  the  introduction  of  a  streamfunction  (u  =  dii'/dy  aud 
v  =  —di'/dx),  the  equations  undergo  a  transformation,  which  scales  the  normal  coordinate  y.  streamfunction 
t,  and  the  external  velocity  ue  with  reference  to  a  constant  velocity  u<>  and  the  local  streamwise  coordinate 
x 


/"o 


/(*,?)  = 


**(*.«)- 


y/U0UX 


4A*iV) 


where  «o  is  taken  as  the  magnitude  of  vector  mean  velocity.  The  concept  of  constant  boundary  layer 
thickness,  attained  by  Falkner-Skan  variables  (with  ue  as  reference  velocity),  has  to  be  abandoned  because 
the  external  velocity  is  unknown  in  interactive  calculations.  However,  provided  that  the  integration  does  not 
start  in  the  immediate  neighbourhood  of  the  stagnation  point,  growth  of  the  boundary  layer  thickness  can 
be  kept  limited.  In  terms  of  these  so-called  semi-transformed  coordinates  equations  .written  as  a  first  order 
system  by  means  of  two  additional  dependent  variables  U  and  V,  and  boundary  conditions  take  the  form 

U*  -  V 

2  dx  V     Ox  Ox' 
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The  discretization  of  the  flow  field  follows  closely  the  above  outlined  procedure  of  the  direct  method,  covering 
the  generation  of  au  orthogonal  grid  and  the  introduction  of  central  differences  and  two-point-averages.  The 
overall  solution  proceeds  in  the  downstream  direction,  accounting  for  downstream  travelling  disturbances 
only.  In  regions  of  backflow  disturbances  propagate  against  the  direction  of  integration  causing  numerical 
instabilities.  On  the  assumption  that  backflow  velocities  are  comparatively  small,  a  stable  integration  can 
be  carried  out  by  adopting  the  Flare  approximation  (Flugge-Lotz  and  Reyhner).  in  which  the  streamwise 
convection  term  udu/dx  is  set  equal  to  zero  in  regions  of  backflow.  With 
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designating  an  "on-off-switch"  of  the  streamwise  convection  term,  the  finite  difference  equations  of  the 
interactive  boundary  layer  problem  become 


&£%Z&.*U^^{^ZL^ 
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Boundary  conditions  are  expressed  in  terms  of  nodal  values  ,too,  whereby  the  evaluation  of  the  integral 
occuring  in  the  interactive  boundary  condition  involves  an  approximation  in  the  fashion  of  the  panel  method 
approach,  leading  to 

U{  m  0,         ft  =  0 

U'j=W},      Wi^h  +  etiiWln-fi) 

where  g,  and  c,,  denote  a  parameter  and  the  diagonal  element  of  the  interaction  matrix,  resulting  from 
a  discrete  approximation  to  the  Hilbert  integral.  Averaging  as  well  as  centering  is  supposed  to  obey  the 
principle  that  the  number  of  generated  terms  approaches  a  minimum.  This  entails  ordinary  differential 
equations,  like  the  y-momentum  equation,  being  centered  about  the  middle  of  the  downstream  face  and 
partial  differential  equations,  like  the  x-momentum  equation,  being  centered  about  the  middle  of  the  box 

A  balance  of  unknowns,  which  occur  here  as  vectors  in  four  components.  {/J,  (7J-,  Vj.Wj}T,  confirms 
the  principal  solvability  of  the  system:  The  J  quadruplets  of  unknowns  match  with  a  total  of  4  J  equations, 
including  2(J-1)  auxiliary  relations.  (J—  1)  x-momentum  and  [J— I)  y-momentum  equations,  each  of  which 
corresponding  to  one  of  the  ( J  —  1)  net  rectangles,  and  4  boundary  conditions.  After  linearizing  this  system 
around  the  values  of  the  preceding  iteration  (iteration  counter  "k— 1"),  respectively  around  the  solution  of 
the  adjacent  upstream  x-station  in  case  of  the  first  iteration,  we  arrive  at  a  linear  system  in  the  Newton 
iterates  6 f*f* ,  SU'f* ,  6VJ-* ,  5W'/K 
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supplemented  by  the  two  components  of  no-slip  condition,  edge  and  interactive  boundary  condition 
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Terms  have  been  grouped  such  that  known  quantities  reside  on  the  right  hand  sides,  while  unknown  quantities 
appear  left  of  the  equal  sign.  The  abbreviated  coefficients  in  the  momentum  equation  are  defined  by 


Mi   =    h~~ +  ~2kT  (/j^2  "  **-*'     lfj 
[*7) j    ■-    gj    +    2A%.     U,-i/a      Jj-i/t)  +  4/y-i 

/.    W.K fi^l/8  rri.K-l    9, 


and  the  right  hand  side  of  the  very  equation  reads 


-1/2 


+       fe.        [{"     )j-i/2       \L/     i;-l/2'y;-l/2  +  V;-l/2   /j-Va    +    Vj-l/2;j-l/2         ^j-1/2  V;-l/2jj 
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Since  the  overall  procedure  involves  a  repetitive  linear  pattern  to  approach  the  solution  of  the  nonlinear 
system,  the  linear  equation  solver  has  to  be  as  fast  as  possible.  Fast  algorithms  take  advantage  of  specific 
matrix  structures,  oue  of  which  is  the  block-tridiagoual  structure,  which  can  be  enforced  in  this  method  by 
the  following  arrangement  of  equations  and  boundary  conditions: 

a.  First  set  of  equations  (index  j  =  1) 

1.  Wall  boundary  condition  prescribing  no  penetration 

2.  Wall  boundary  condition  prescribing  zero  tangential  velocity 

3.  Second  auxiliary  relation  (linking  U'  and  V")  of  the  first  box 

4.  Trivial  y-momentum  equation  of  the  first  box  first  box 

b.  Intermediate  sets  of  equations  (index  2  <  j  <  J  -I) 

1.  First  auxiliary  relation  (linking  /'  and  U)  of  the  (j-'-l)-st  box 

2.  Momentum  equation  of  the  [j—  l)-st  box 

3.  Secoud  auxiliary  relation  (linking  U'  and  V")  of  the 

4.  Trivial  y-momentum  equation  of  the  y-th  box  y-th  box 

c.  Last  set  of  equations  ( index  j  =  J) 

1.  First  auxiliary  relation  (linking  /'  and  U)  of  the  last  box 

2.  Momentum  equation  of  the  last  box 

3.  Interactive  boundary  condition 

4.  Edge  boundary  condition. 
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Following  these  instructions  equations  and  boundary  conditions  can  be  summarized  in  matrix-vector  form, 
offering  the  desired  structure 
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where  {#*}  =  {i^^iYJ^.tWf^  and  {rf }  =  [trJfA'df'.l'M'Au)?}*  denote  four  di- 
mensional  vectors  of  the  unknown  Newton  iterates  and  the  known  right  hand  sides,  respectively.  The  blocks 
in  the  three  diagonals  of  the  above  matrix  have  the  dimension  4x4  and  can  be  obtained  from 
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The  components  of  the  right  hand  side  vectors  can  be  calculated  from  the  following  formulae 
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The  numerical  solution  of  the  above  system  can  again  be  achieved  by  Kellers  block  elimination  method, 
which  works  very  much  like  a  Gauss's  algorithm,  but  .firstly,  matrices  are  eliminated  instead  of  scalars.  and. 
secondly,  quite  a  few  manipulations  can  be  saved  because  of  sparse  occupation. 
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3.3  Interaction  model 

The  coupling  of  the  boundary  layer  and  external  inviscid  flow  is  referred  to  as  interaction  modeL  Basically, 
the  interaction  consists  of  thickening  the  effective  blade  shape  by  viscous  displacement,  which  in  turn  causes 
changes  in  the  surface  pressure.  The  exchange  of  information  between  viscous  and  inviscid  regions  takes 
place  via  the  boundary  conditions  only,  which  is  the  specification  of  either  an  impermeable  displacement 
surface  or  a  nonzero  wall  transpiration  at  the  original  surface  in  case  of  inviscid  flow,  respectively  the 
prescription  of  either  pressure,  or  displacement  thickness,  or  a  linear  combination  of  both  in  case  of  viscous 
flow.  The  interaction  is  called  weak  if  viscous  effects  on  pressure  remain  smalL  while  situations  where 
viscous  disturbances  to  the  inviscid  flow  field  are  substantial  demand  the  application  of  strong  simultaneous 
interaction.  Interaction  models  fall  into  four  types,  the  first  three  of  which  provide  loose  coupling  of  viscous 
and  inviscid  regions: 

1.  The  direct  interaction  method  combines  a  direct  inviscid  and  a  direct  viscous  flow  solver.  This  classi- 
cal approach  achieves  a  solution  by  iterative  matching  of  boundary  layer  calculations  with  prescribed 
pressure  and  inviscid  calculations  with  prescibed  displacement  thickness.  The  alternate  treatment  of 
pressure  and  displacement  thickness  leads  to  a  local  breakdown  of  the  procedure  when  slight  changes 
in  pressure  entail  significant  changes  in  displacement  thickness. 

2.  In  inverse  methods  the  roles  of  displacement  thickness  and  pressure  are  interchanged.  Hence  the  in- 
viscid flow  equations  determine  the  displacement  thickness  distribution,  which  is  imposed  as  boundary 
condition  on  the  viscous  flow  calculation.  The  results  of  which  is  the  pessure,  which  closes  the  cycle 
by  being  input  to  the  inviscid  flow  computation.  The  hierarchical  manner  of  solving  the  complete  flow 
problem  excludes  this  model  just  as  the  previous  one  from  handling  strongly  interacting  regions. 

3.  The  semi-inverse  interaction  method  is  composed  of  a  direct  inviscid  and  an  inverse  viscous  flow  solver. 
Both  parts  of  the  scheme  process  the  same  input,  ie.  displacement  thickness,  and  generate  the  same 
output,  ie.  pressure.  After  each  cycle  an  updated  displacement  thickness  is  relaxed  based  on  the 
deviation  of  the  two  pressure  distributions,  which  should  coincide  upon  convergence. 

The  existence  of  a  definite  hierarchy  between  boundary  layer  and  outer  flow  in  the  above  mentioned  methods 
gives  rise  to  a  limited  rate  of  feedback  between  both  regions.  Strong  interaction  models  incorporate  the  outer 
flow  somehow  in  the  boundary  layer  calculations,  foi  example  by  the  following  interaction  law: 

4.  Simultaneous  (or  strong)  interaction  methods  solve  the  boundary  layer  equations  subject  to  an  interac- 
tion law  as  outer  boundary  condition,  which  retrieves  the  elliptic  character  of  the  outer  flow.  No  definite 
assignment  of  displacement  thickness  and  pressure  can  be  made  to  viscous  and  inviscid  region.  Rather, 
both  quantities  are  treated  as  unknowns,  related  by  the  interaction  law.  The  procedure  emphasizes 
simultaneous  solution  for  both  displacement  thickness  and  pressure. 

In  the  following  the  currently  used  interaction  law  will  be  formulated  in  terms  of  displacement  thickness 
and  external  velocity,  the  latter  being  related  to  pressure  by  Bernoulli's  equation.  Both  inviscid  external 
velocity,  resulting  from  a  potential  flow  calculation  around  the  original  blade,  ueI,  and  a  correction  term  due 
to  viscous  displacement,  ues,  are  assumed  to  contribute  to  the  total  external  velocity 

ue(x)  =  ueI{x)  +  ue6(jr) 

The  contribution  due  to  viscous  effects  can  be  evaluated  by  means  of  the  blowing  velocity  concept,  which 
has  its  origin  in  a  paper  by  Lighthill  [16].  He  proved  that  the  effect  of  boundary  layers  on  the  outer  flow  can 
be  represented  by  a  surface  distribution  of  sources.  Nonzero  normal  velocities  at  the  surface  of  blades,  being 
induced  by  this  very  source  distribution,  are  intended  to  displace  streamlines  from  the  surface  such  that  the 
virtual  displacement  surface  becomes  a  streamline  (see  figure  4  on  page  12) 

dfjx)       vjx.6*) 

dx  Ue(*) 

where  v(x. 6*)  denotes  the  vertical  velocity  component  on  the  displacement  surface.  Since  it  is  sufficient 
to  approximate  the  correction  term  u«$,  further  deductions  make  use  of  the  thin  airfoil  approach.  Firstly, 
streamwise  upper  and  lower  surfaces  of  the  blade  will  be  considered  flat  plates  in  this  context,  implying 
that  the  blowing  velocity  equals  half  the  local  source  strength,  and  , secondly,  the  displacement  thickness  is 
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supposed  to  be  so  small  that  the  horizontal  velocity  components  of  the  inviscid  flow  do  not  vary  across  the 
boundary  layer 


lf-«(*,o).v{xX)-l  £*"5<"''*> 


Unlike  vertical  components  horizontal  velocities  are  affected  by  all  sources  distributed  along  the  supposedly 
flat  surface.  Eventually  the  viscous  correction  term  is  given  by  the  well  known  thin  airfoil  integral,  which  is 
referred  to  here  as  Hilbert  integral 
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As  indicated  by  the  integration  limits  interaction  is  restricted  to  a  finite  region  *&<*<*«  on  either 
streamwise  upper  or  lower  surface.  The  numerical  implementation  of  the  interaction  law  requires  some 
discrete  approximation  of  the  above  thin  airfoil  integral.  Adopting  the  panel  method  approach,  which 
concerns  here  a  piecewise  approximation  of  the  continuous  blowing  velocity  d{ueS*)/dx  to  allow  piecewise 
analytical  integration,  the  integral  can  be  written  as  a  finite  series 
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with  [cue]  denoting  a  matrix  of  interaction  coefficients.  Recalling  that  the  boundary  layer  calculation  at  a 

streamwise  position  involves  iterations,  the  inviscid  contribution  can  be  included  in  the  total  external  velocity 

of  the  previous  Newton  iteration  (identified  by  "k—  1"),  leading  to  a  generalized  version  of  the  interaction 

law 

K 

Since  displacement  thickness  does  not  belong  to  the  dependent  variables.  (<!>*)'•*,  the  displacement  thickness 
of  the  streamwise  position  currently  being  solved,  must  be  expressed  in  terms  of  dependent  variables  to  make 
allowance  for  its  unknown  status.  With  (6*)''K  being  replaced  by  yj  —  ^jK/[ue)%'K  aU(l  after  separating 
known  and  unknown  terms,  the  interaction  law  takes  its  almost  final  form 

(««)'•*  (1 -<?««/ )  +  <**/*-# 

with 

g?  =  [«<Y-K-X  +   t  <**[(«€*•)*•*  -  (M')*'*""1]  -  <*<(«.«/  -  op""1 

Implementing  this  relation  necessitates  a  known  right  hand  side,  which  can  be  evaluated  only  in  an  ap- 
proximate fashion,  because  the  terms  [ue6*)k'K  are  not  known  yet  downstream  of  the  current  x-location. 
To  ensure  interaction  also  over  this  region,  these  terms  are  taken  from  the  previous  iteration  updated  by 
some  relaxation  formula.  With  e;fc  ?=  e/fc  y/vx/uo  and  gf  ss  gf/uo  designating  dimensionless  interaction 
coefficients  and  right  hand  side,  respectively,  the  actually  coded  interaction  law  can  be  written  in  terms  of 
semi-transformed  coordinates 

This  equation,  being  used  as  an  outer  boundary  condition  in  the  viscous  flow  computation,  relates  the 
unknown  external  velocity  and  the  unknown  displacement  thickness  of  the  i-th  boundary  layer  station, 
whereby  displacement  thickness  has  been  expressed  by  streamfunction  and  external  velocity.  Because  of  the 
elliptic  character  of  the  outer  flow,  which  has  been  incorporated  in  the  boundary  layer  via  the  interaction 
law.  solution  requires  a  global  iteration,  consisting  of  several  viscous  sweeps  to  be  performed  over  both  the 
streamwise  upper  and  lower  surface. 
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Note  on  the  calculation  of  interaction  coefficients:  Let  the  continuous  function  "external  velocity  times 
displacement  thickness",  denoted  by  D,  be  discretized  at  a  finite  number  of  streamwise  positions.  The  thin 
airfoil  integral  then  can  be  approximated  by  a  finite  series  of  weighted  "DV  at  the  very  locations 
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1   [**dD     dt  £ 


The  weights  are  the  interaction  coefficients  c,*,  with  the  first  subscript  indicating  the  streamwise  position, 
where  the  correction  term  ttet  is  to  be  evaluated,  and  the  second  indicating  the  location,  whose  effect  is 
accounted  for.  With  a  properly  interpolated  D-function,  integration  can  be  carried  out  analytically  piece 
by  piece.  Provided  the  point  under  consideration  does  not  fall  within  the  limits  of  integration.  D  will  be 
interpolated  linearly.  Piecewise  linear  functions  can  be  built  up  by  overlapping  triangular  distributions, 
integration  over  which  yields  the  coefficients  whose  k  3=  f,  k  ft  i—1,  k  £  i+l 
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A  linearly  interpolated  D-function  would  lead  to  singular  integrals  for  the  coefficients  k  =  1,  k  =  i— 1,  and 
k  SB  i  +  l.  Therefore  D  will  be  approximated  by  a  polynomial  of  degree  2  in  the  interval  a_i  <  £  <  Xi+i  . 
Splitting  again  into  overlapping  distributions,  which  this  time  are  parabolic,  and  applying  Cauchy's  principal 
value  technique  permits  integration  of  singular  integrands: 
The  coefficient  at  the  middle  of  the  inducing  source  distribution  is  given  by 
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The  coefficient  at  the  right  of  the  inducing  source  distribution  is  given  by 
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The  coefficient  at  the  left  of  the  inducing  source  distribution  is  given  by 
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for    j-,_i  <  i  <  xi+l 

for    xi+l  <  i  <  x,+2 
•r.-A+i 
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As  indicated  above,  the  overall  solution  is  approached  in  an  iterative  process,  in  which  alternately 
viscous  and  inviscid  flow  equations  are  being  solved  (see  also  figure  1  on  page  4): 

1.  Calculate  the  external  velocity  distribution  ue/  in  an  inviscid  flow  field  by  means  of  the  panel  method. 
To  account  for  the  blade-thickening  due  to  viscous  displacement,  specify  a  nonzero  normal  velocity 
( blowing  velocity)  at  the  blade's  surface. 

2.  March  through  the  boundary  layers  of  both  streamwise  upper  and  lower  surface  using  (over  the  main 
part  of  the  blade)  the  interaction  law  as  outer  boundary  condition. 

3.  Check  convergence  and  quit  if  the  criterion  is  satisfied. 

4.  If  the  convergence  criterion  is  not  met.  prepare  for  another  cycle.  Update  the  product-term  "external 
velocity  times  displacement  thickness**  on  the  basis  of  the  deviation  between  inviscid  and  viscous  external 
velocity  distribution 


(«/)^s(MfA[l  +  «(ffja-l)' 


where  A  denotes  the  counter  of  global  iterations.  Further,  compute  the  blowing  velocity  distribution, 
which  serves  as  boundary  condition  for  the  inviscid  flow  solution 
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and  proceed  with  the  first  step. 


3.4  Turbulence  modelling 

The  presence  of  additional  unknown  shear  stresses  in  turbulent  flows  requires  modelling  assumptions  to 
balance  the  number  of  unknowns  and  equations.  Eddy  viscosity  models,  one  of  which  is  used  in  the  present 
method,  relate  turbulent  shear  stresses  to  mean  flow  quantities  on  an  empirical  basis.  They  draw  their 
versatility  from  the  convenience  of  maintaining  the  same  approach  and  numerical  formulation  for  both 
laminar  and  turbulent  flows.  The  currently  employed  model  features  a  zero  equation  closure,  which  presumes 
that  the  length  scale  can  be  specified  by  an  algebraic  equation.  The  turbulent  kinematic  viscosity  is  defined 
by  two  separate  formulae,  one  for  the  inner  region  being  based  on  Van  Driest's  approach  [21],  and  the  other 
for  the  outer  region  being  based  on  a  velocity  defect  approach 
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Continuity  of  the  turbulent  kinematic  viscosities  is  established  by  defining  yc  as  the  distance  from  the  wall 
where  expressions  for  inner  and  outer  region  do  agree.  Since  transition  is  not  an  instantaneous  process, 
an  intermediate  status  of  flow,  called  transitional  flow,  is  assumed  between  the  laminar  and  fully  turbulent 
regions.  This  region  is  taken  into  account  by  introduction  of  an  intermittency  factor,  which  smears  out  the 
step-shaped  change  from  kinematic  to  eddy  viscosity 


where  RXfr  denotes  the  Reynolds  number  based  on  external  velocity  and  streamwise  location  xfr  at  onset 
of  transition.  A  value  of  1200  was  originally  assigued  to  the  empirical  constant  G-!(r.  However,  numerical 
experiments  seem  to  indicate  that  low  Reynolds  number  flows  can  be  better  modelled  by  values  below  1200. 
The  parameter  a  in  the  outer  region  formula  is  obtained  from 
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where  the  nondimensional  factor  F  denotes  the  ratio  of  total  turbulence  energy  production  to  the  shear- 
stress-related  turbulence  energy  production,  evaluated  at  the  location  of  maximum  shear  stress 
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The  ratio  of  the  time-averaged  quantities  above  can  be  approximated  by  a  function  of  Rt  =  r«,/(  — uV), 
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Since  this  turbulence  model  is  not  validated  for  separated  flow,  eddy  viscosities  in  regions  of  backflow 
correspond  to  those  of  the  adjacent  upstream  station  at  the  same  ij -coordinate. 

Transition,  if  not  available  from  other  sources,  currently  is  predicted  by  an  empirical  data  correlation 
expressed  in  terms  of  the  Reynolds  numbers  based  on  momentum  thickness  and  streamwise  coordinate  at 
onset  of  transition 

n*tr 

Should  laminar  transition  occur  upstream  of  the  computed  point  of  transition,  then  transition  is  redefined 
at  the  point  of  laminar  separation. 

4.  SAMPLE  CALCULATIONS 

This  chapter  will  focus  ou  sample  calculations  reporting  on  the  influence  of  turbulence  model  and  loci  of 
transition,  as  well  as  ou  the  sensitivity  of  calculations  to  the  choice  of  the  relaxation  parameter.  Further 
calculations,  referring  to  the  variation  of  air  inlet  angle  and  Reynolds  number,  are  intended  to  demonstrate 
the  code's  abilities  to  predict  the  phenomena  of  low  Reynolds  number  and  high  angle  of  attack  flow  regimes. 
No  comparisons  with  experimental  data  will  be  presented  in  this  report.  However,  the  severe  deviation  of 
early  calculations  from  experimental  results,  together  with  numerical  breakdowns  and  convergence  problems, 
revealed  the  need  for  a  modification  of  the  turbulence  model  in  the  low  Reynolds  number  regime.  Accord- 
ingly, the  transition  length  has  been  reduced,  resulting  in  fewer  numerical  troubles  and  improved  agreement 
of  computational  and  experimental  results.  Since  the  development  of  the  code  is  not  completed  yet  and 
further  improvements  are  to  be  expected,  the  final  evaluation,  including  the  presentation  of  comparisons 
with  experiments,  is  still  to  come.  All  the  results  describe  flows  past  cascades,  whose  blades  are  staggered  by 
an  angle  of  45  degrees  and  whose  solidity  amounts  to  0.78.  Cross  sections  of  the  NACA-65  series  compressor 
blade  family  were  selected  in  consideration  of  available  experimental  data.  In  the  following,  results  of  test 
case  computations  and  a  brief  analysis  of  them  will  be  given. 

Figures  5.  6.  7.  and  S  show  the  convergence  histories  as  a  function  of  the  relaxation  parameter  u/.  which 
has  a  major  influence  on  the  number  of  iterations  required  to  obtain  converged  results.  Relaxation  is  meant 
to  ensure  or  to  accelerate  convergence.  The  lift  coefficient  based  on  inlet  velocity  is  taken  as  convergence 
indicator,  although  other  choices  are  equally  possible.  Figure  5  presents  the  convergence  behaviour  of  a 
high  Reynolds  number  flow  at  zero  incidence.  The  geometrical  input  consists  of  51  blade  coordinates,  taken 
from  the  original  NACA  publication.  The  code  works  with  recomputed  coordinates,  which  are  cosine- 
distributed  along  the  chord  length.  Subplot  (a)  refers  to  51  of  these  recomputed  coordinates,  subplot  (b) 
to  97.  The  upper  plot  of  figure  5  displays  a  desirable  convergence  behaviour,  wherein  all  curves  approach 
the  same  lift  coefficient  after  a  finite  number  of  iteration  cycles.  Maximum  acceleration  is  achieved  by  a 
relaxation  parameter  of  2.0  or  slightly  higher.  Smaller  relaxation  parameters  delay  convergence,  higher 
relaxation  parameters  can  lead  to  numerical  breakdowns.  The  lower  plot  of  figure  5  indicates  a  numerical 
instability,  which  seems  to  be  related  to  the  panel  code.  There  is  no  evidence  that  these  oscillations  are  due 
to  unsteady  loci  of  transition.  The  oscillations  of  the  lift  coefficients  are  stable  in  amplitude  and  the  mean 
of  all  oscillations  tends  towards  a  single  value.  Since  the  amplitude  of  the  oscillations  amounts  only  to  half 
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a  percent  of  the  lift  coefficient,  results  might  even  be  regarded  as  converged.  Probably  not  the  panel  code 
itself,  but  the  procedure  how  coordinates  are  recomputed  is  responsible  for  the  instabilities.  The  currently 
used  interpolation  functions  are  -  unlike  spline  functions  -  continuous  in  values  of  functions  only,  not  in 
derivatives.  If  the  number  of  recomputed  coordinates  strongly  exceeds  the  number  of  input  coordinates, 
then  the  used  interpolation  algorithm  might  give  rise  to  irregular  shapes  in  the  leading  and  trailing  edge 
areas,  where  points  are  packed  because  of  the  cosine-redistribution. 

The  same  cascade  configuration  at  the  same  incidence,  but  at  a  much  smaller  Reynolds  number  leads  to 
numerical  divergence  and  numerical  termination  of  the  calculation,  as  can  be  seen  from  figure  6.  Subplots 
(a)  and  (b)  refer  again  to  51  and  97  recomputed  coordinates.  Calculations  have  been  performed  up  to 
40  iteration  cycles,  of  which  only  the  first  20  are  shown.  Some  of  the  computations  break  down  between 
cycle  20  and  40,  the  remaining  computations  keep  changing,  but  do  not  get  terminated  within  40  cycles. 
Furthermore,  computations  for  different  relaxation  parameters  do  not  approach  a  common  lift  coefficient. 
Numerical  breakdowns  are  either  due  to  a  failure  of  the  Newton  iterations  at  a  particular  streamwise  position 
to  converge,  or  due  to  negative  external  velocities  as  an  artificial  result  of  extensively  separated  boundary 
layers.  These  difficulties  are  related  to  unrealistically  large  regions  of  flow  separation,  which  are  caused  by 
the  turbulence  modeL  The  original  CS-modeL  being  based  on  high  Reynolds  number  data,  overestimates 
the  length  of  transitional  flow  for  low  Reynolds  numbers,  which  in  conjunction  with  bubble  transition  can 
trigger  dramatic  deviations  from  realistic  flow  results. 

Transitional  flow  in  the  CS-model  is  accounted  for  by  an  intermittency  factor  itr 


Itr  -  1  -  <?*p 


This  factor,  being  a  transient  function  with  values  0  at  onset  of  transition  and  1  at  fully  turbulent  flow, 
avoids  an  abrupt  transition  from  laminar  to  turbulent  flow  by  smoothing  out  the  step-shaped  change  from 
kinematic  to  eddy  viscosity.  The  length  of  transitional  flow  is  mainly  determined  by  the  empirical  constant 
G^(r,  which  was  set  1200  in  the  original  CS-model.  This  choice  makes  the  transitional  region  of  low  Reynolds 
number  flows  spreading  too  far  into  the  turbulent  region.  Consequently,  low  Reynolds  number  flows  with 
a  G'7fr  of  1200  are  more  likely  to  separate  and  less  likely  to  reattach  in  the  transitional  regiou.  Numerical 
experiments  (see  also  figures  9  to  15)  imply  that  low  Reynolds  number  flows  require  smaller  constants  G-!fr. 
Figure  7  shows  the  results  of  the  same  low  Reynolds  number  flows  as  in  figure  6.  using  a  slightly  modified 
turbuienec  model,  but  still  computing  the  point  of  transition.  Subplot  (a)  corresponds  to  a  G~.fr  of  120. 
subplot  (b)  to  a  G-,fr  of  12.  respectively.  Calculations  seem  to  converge  now.  however  at  retarded  rates  of 
convergence  when  compared  with  the  high  Reynolds  number  case.  Most  of  the  computations  end  up  with 
oscillations  of  the  lift  coefficient,  which  originate  here  from  moving  transition  points.  Fixing  transition  on 
both  upper  and  lower  surface  turns  the  oscillating  lift  histories  of  figure  7  into  smooth  curves,  which  are 
given  in  figure  8.  Transition  on  the  uppper  surface  is  assumed  at  midchord.  transition  on  the  lower  at  70 
percent  of  chord.  Fixing  transition  supresses  a  degree  of  freedom,  because  of  which  only  carefully  selected 
locations  might  work  out.  The  above  chosen  loci  reflect  the  most  upstream  positions,  where  trausition  ever 
has  been  experienced  within  20  cycles.  Recapitulating,  reasonable  results  for  a  low  Reynolds  number  flow 
could  be  accomplished  in  the  investigated  case  by  the  following  provisions: 

1.  The  extension  of  transitional  flow  has  to  be  narrowed  by  decreasing  an  empirical  constant  in  the  CS- 
model. 

2.  The  loci  of  transition  have  to  be  fixed  at  reasonable  locations. 

The  following  6  plots  offer  some  insight  into  the  convergence  histories  of  the  low  Reynolds  number  flows, 
whose  lift  histories  are  given  in  figure  C.a  and  8.a.  respectively.  Figures  9.  11.  and  13  refer  to  the  unsuccessful 
calculation,  using  the  original  CS-model.  computed  transition,  and  a  relaxation  parameter  of  2.  Figures  10. 
12.  and  14  display  the  results  of  an  equivalent  calculation,  but  now  using  the  modified  turbulence  model  with 
a  Gftr  of  120.  fixed  transition  at  midchord  on  the  upper  and  70  percent  of  chord  on  the  lower  surface,  and  a 
relaxation  parameter  of  2.  The  histories  of  pressure  coefficient,  local  skin  friction  coefficient,  and  displacement 
thickness  are  presented  for  the  first  10  cycles.  Figure  9  and  10  give  the  histories  of  pressure  distributions, 
showing  the  developments  from  purely  inviscid  results  (indicated  by  iteration  number  0)  to  those  altered 
increasingly  by  viscous  displacement  effects.   Pressure  distributions  are  rather  insensitive  compared  to  the 
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more  dramatic  changes  of  the  lift  coefficients  during  the  iteration  process.  With  the  exception  of  the  trailing 
edge,  pressure  distributions  keep  looking  the  same,  only  flow  separation  leaves  some  marks  -  a  sort  of  flat  - 
on  the  pressure  distributions.  The  longer  flats  in  figure  9  suggest  a  more  extended  flow  separation  than  in 
figure  10. 

Figures  11  and  12 ,  which  show  the  developments  of  local  skin  friction  distributions,  prove  this  assertion. 
Two  flow  phenomena  can  be  read  from  skin  friction  distributions:  Flow  separation  occurs  at  the  point  of  zero 
skin  friction  and  transition  is  accompanied  by  a  steep  increase  in  skin  friction.  A  closeup  of  the  numerically 
unstable  computation  of  figure  11  indicates  that  laminar  flow  separation  takes  place  a  few  percent  upstream 
of  midchord  on  the  upper  surface,  but  turbulent  reattachment,  which  closes  the  bubble  usually  not  too 
far  downstream  of  the  point  of  separation,  seems  to  be  omitted.  At  the  streamwise  position,  where  fully 
turbulent  flow  is  established,  backflow  dominates  so  that  the  boundary  layer  is  not  able  to  return  to  the  wall. 
The  lower  surface  exhibits  laminar  separation  too,  but  there  the  flow  does  form  a  bubble.  These  results  are 
considered  unsatisfactory,  because  one  would  not  expect  a  severe  trailing  edge  separation,  extending  over 
half  the  chordlength,  in  cascades  with  only  slightly  cambered  blade  sections  at  zero  incidence,  even  when 
exposed  to  low  Reynolds  numbers.  The  results  of  figure  12.  which  were  obtained  by  using  the  modified 
turbulence  model  and  fixed  transition,  fit  much  better  both  our  expectations  and  experiments.  Forced 
transition  makes  the  results  in  some  way  less  realistic,  but  forced  transition  is  necessary  to  avoid  oscillations. 
Transition  is  followed  on  both  surfaces  by  a  steep  increase  in  skin  friction,  being  caused  by  the  fuller  velocity 
profiles  of  turbulent  boundary  layers.  The  locations  of  transition,  which  had  been  chosen  specifically  to 
circumvent  numerical  instabilities,  seem  to  suppress  laminar  separation  and  bubble  transition.  The  flow 
remains  completely  attached  on  the  lower  surface,  while  it  separates  on  the  upper  surface  close  to  the 
trailing  edge.  This  turbulent  trailing  edge  separation  can  be  noticed  for  the  first  time  in  cycle  3.  and  a  few 
cycles  afterwards  a  steady  separation  over  6  percent  of  chordlength  has  been  established. 

Displacement  thickness  measures  the  deviation  of  viscous  flow  from  inviscid  flow.  i.e.  it  tells  how  far 
streamlines  are  shifted  from  the  surface  by  viscous  effects.  Uneven  displacement  on  upper  and  lower  surface 
causes  a  lift  loss  compared  to  the  inviscid  flow  case.  Figure  13  and  14.  which  present  the  developments 
of  displacement  thickness  distributions,  depict  the  differences  of  displacement  thickness  on  the  upper  and 
lower  surface,  which  tend  to  uncamber  the  airfoil.  Displacement  thickness  at  the  trailing  edge  -  like  the 
lift  coefficient  -  proves  to  be  an  excellent  convergence  indicator.  The  breakdown  of  calculations  in  figure 
13  is  preceded  by  an  impressive  growth  of  displacement  thickness  at  the  trailing  edge  on  the  upper  surface. 
Constant  values  since  cycle  8  signal  convergence  of  that  case  shown  in  figure  14.  Both  figures  have  marks, 
left  by  transition  aud  separation.  Transition  is  usually  accompanied  by  slight  drawback  in  the  steady  growth 
of  displacement  thickness,  which  can  be  seen  clearly  on  both  surfaces.  Separation  dramatically  increases 
the  growth  of  the  displacement  thickness,  which  is  well  demonstrated  at  the  trailing  edge  on  both  upper 
surfaces. 

Figure  15  compares  the  results  for  two  different  airfoil  sections  of  the  NACA  G5-series  compressor  Made 
family,  both  of  which  are  arranged  in  a  staggered  cascade  of  relatively  low  solidity.  The  variations  of  lift 
coefficient  based  on  inlet  velocity  and  turning  angle  as  a  function  of  air  inlet  angle  are  given  for  both  inviscid 
and  converged  viscous  flow  computations.  The  blade  sections  are  the  10  percent  thick  NACA  C5-410  and 
NACA  C5-(  12)10.  whose  camber  is  indicated  by  the  first  number  after  the  dash,  denoting  the  design  lift 
coefficient  in  tenths  for  the  isolated  airfoil.  The  high  Reynolds  number  did  allow  calculations  to  he  carried 
out  with  the  original  CS-turbuleuce  model  and  having  the  points  of  transition  computed.  In  spite  of  the 
high  Reynolds  number  inviscid  and  viscous  results  are  clearly  different  due  to  viscous  displacement  effects. 
which  are  more  pronounced  in  case  of  the  more  cambered  blade  section.  The  beginning  breakaway  of  the 
viscous  curves  from  their  inviscid  counterparts,  observed  in  the  left  subplot,  suggests  ouset  of  stall  at  an  air 
inlet  angle  of  62  degrees  for  the  investigated  cascade.  Calculations  for  air  inlet  angles  at  and  beyond  stall 
could  not  be  performed  successfully.  Therefore  the  plots  were  terminated  before  approaching  angles  of  deep 
stall. 

The  following  three  plots  reveal  the  dependence  of  local  flow  quantities  on  air  inlet  angle  for  that  cascade 
flow,  whose  lift  variations  are  given  in  the  left  plot  of  figure  15.  The  variations  of  pressure,  skin  friction,  and 
displacement  thickness  -  all  in  dimensiouless  form  -  are  shown  in  figure  16.  17  and  18.  The  distributions  of 
pressure  coefficients,  based  on  the  dynamic  pressure  of  inlet  velocity,  provide  information  about  the  blade 
loadings.  Pressure  peaks  result  from  highly  accelerated  flow  around  the  leading  edge,  when  the  stagnation 
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point  is  relatively  far  away  from  the  geometric  leading  edge.  Only  those  pressure  distributions,  which 
correspond  to  air  inlet  angles  around  the  design  lift  coefficient,  can  avoid  velocity-peaks  at  the  leading  edge. 
The  flattening  of  pressure  distributions  in  the  trailing  edge  area  indicates  separation  on  the  upper  surface 
for  the  highest  incidence  angles.  None  of  the  flows  exhibits  a  laminar  separation  on  either  upper  or  lower 
surface,  which  anyway  would  not  be  expected  for  such  a  high  Reynolds  number.  At  air  inlet  angles  of  60  and 
62  degrees  turbulent  trailing  edge  separation  (without  reattachment)  does  occur,  extending  over  4  and  20 
percent  of  chordlength,  respectively.  Transition  from  laminar  to  turbulent  flow  regime  is  accompanied  by  a 
step-shaped  increase  of  skin  friction,  which  permits  the  identification  of  transition  in  the  skin  friction  plots. 
With  increasing  incidence  transition  is  shifted  on  the  upper  surface  towards  the  leading  edge,  while  the  lower 
surface  experiences  a  reverse  move  of  transition.  The  development  of  displacement  thickness  distributions 
gives  a  good  indication  how  much  flow  is  disturbed  by  viscous  effects.  Generally,  the  perturbation  increases 
with  higher  incideuce.  However,  for  positive  incidence  the  lower  surface  shows  no  significant  changes  in  the 
distributions  of  displacement  thickness.  It  is  the  difference  between  the  displacement  thickness  on  upper  and 
lower  surface,  that  entails  the  loss  of  lift.  Since  the  difference  is  more  pronounced  at  high  incidence,  those 
flows  depart  more  strongly  from  the  inviscid  ideal 

Finally,  the  last  four  plots  are  intendend  to  promote  the  basic  understanding  of  viscous  flows  past 
cascades.  Velocity  profiles  along  pressure  and  suction  surface  are  shown  for  four  representative  flow  regimes. 
The  vertical  scales  of  the  velocity  profiles  had  been  enlarged  to  achieve  a  good  resolution,  allowing  the 
illustration  of  marginal  backflow.  Figure  19  and  20  investigate  the  effect  of  the  Reynolds  number  on  viscous 
flow  results.  Flows  at  Reynolds  numbers  of  6  million  and  245,000  were  chosen  to  represent  the  high  and  low 
Reynolds  number  regimes  at  an  air  inlet  angle,  which  corresponds  roughly  to  the  design  lift  coefficient.  Both 
cases  had  been  examined  and  depicted  in  terms  of  other  flow  quantities  on  previous  occasions  in  this  report: 
The  lift  histories  of  the  high  Reynolds  number  flow  in  figure  5. a,  the  lift  histories  and  the  developments  of 
pressure,  skin  friction,  and  displacement  thickness  for  the  low  Reynolds  number  flow  in  figures  8. a.  10, 12.  and 
14.  The  more  obvious  effects  as  the  Reynolds  number  decreases  are  a  more  pronounced  growth  of  boundary 
layers  and  a  shift  of  transition  points  towards  the  trailing  edge.  The  locations  of  transition  in  figure  20 
have  been  specified  (to  avoid  numerical  instabilities);  however,  they  do  reflect  the  most  upstream  positions, 
where  transition  ever  occured  in  the  computation.  The  absence  of  a  somehow  expected  transition  bubble 
in  figure  20  is  clue  to  that  forced  transition,  which,  if  prescribed  upstream  of  the  laminar  separation  point. 
can  eliminate  the  bubble  entirely  or  reduce  its  extent.  Furthermore,  these  figures  demonstrate  the  different 
shapes  of  laminar,  turbulent,  and  separating  velocity  profiles.  The  rapid  flow  acceleration  immediately  at 
the  surface,  as  is  typical  for  turbulent  profiles,  contrasts  the  gradual  flow  acceleration  of  lamiuar  profiles. 

The  influence  of  the  air  inlet  angle  on  boundary  layer  results  is  investigated  in  figure  21  and  22.  Air 
inlet  angles  of  43  and  58  degrees  have  been  selected  to  show  flow  regimes  of  cascades  at  negative  and  positive 
incidence.  As  the  air  inlet  angle  increases,  velocity  profiles  ou  upper  and  lower  surface  differ  more  and 
more,  resulting  in  strongly  asymmetric  displacement  on  the  upper  and  lower  surface.  Transition  on  upper 
surface  moves  with  increasing  air  inlet  angle  towards  the  leading  edge,  while  on  the  lower  surface  the  laminar 
flow  region  spreads  with  increasing  air  inlet  angle,  eventually  leading  to  completely  laminar  flow,  as  can  be 
seen  from  figure  22.  Low  Reynolds  numbers  and  small  or  negative  incidence  angles  can  give  rise  to  bubble 
transition,  which  is  illustrated  in  figure  21.  Flow  ou  upper  surface  remains  laminar  until  midchord.  separates 
then,  becomes  transitional,  and  reattaches  as  a  turbulent  boundary  layer.  Strongly  decelerating  flows  can 
even  cause  the  separation  of  turbulent  boundary  layers.  The  increasing  region  of  backflow  on  the  upper 
surface  at  the  trailing  edge  is  responsible  for  the  loss  of  lift,  as  the  angle  of  attack  approaches  stall.  All  four 
velocity  profile  plots  give  a  good  indication,  how  boundary  layers  grow,  how  they  develop  regions  of  reversed 
flow,  and  also  how  the  external  velocitv  changes  over  the  blade. 
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5.  CONCLUSIONS 

This  report  presents  a  numerical  method  to  analyse  viscous  incompressible  flows  past  2-d  cascades.  Cebeci 
and  collaborators  developed  a  viscous-inviscid  interaction  method,  which  accomplishes  a  complete  flow  so- 
lution by  treating  alternately  inviscid  and  viscous  flow  regimes.  In  viscid  flow  equations  are  solved  by  means 
of  a  standard  panel  method,  -viscous  flow  equations  are  approximated  by  boundary  layer  equations,  whose 
solution  is  determined  by  a  finite  difference  scheme.  The  essence  of  the  method  is  an  interactive  boundary 
condition,  which  forces  viscous  and  inviscid  flow  to  merge  smoothly  at  the  boundary  layer  edge.  It  is  the 
numerical  implementation  of  this  condition,  which,  by  prescribing  a  linear  combination  of  external  velocity 
and  displacement  thickness,  avoids  the  difficulties  incurred  when  either  direct  or  inverse  boundary  layer 
methods  are  employed.  The  current  method  is — like  Navier-Stokes-solvers — capable  of  computing  viscous 
flows  including  flow  separation,  but  in  contrast  to  Navier-Stokes-solvers.  at  significantly  lower  computing 
costs. 

The  potential  and  capabilities  of  the  method  have  been  demonstrated  for  flows  past  cascades  of  moder- 
ately cambered  blade  sections  over  a  range  of  Reynolds  numbers  and  air  inlet  angles.  The  investigated  flows 
involved  both  laminar  separation  bubbles  and  turbulent  separation  without  reattachment.  The  procedure 
seems  to  work  reliably  for  high  Reynolds  numbers  and  small  incidence  angle,  but  encounters  difficulties  in 
cases  of  low  Reynolds  numbers  and  highly  loaded  blades.  Numerical  experiments  indicated  the  need  for  a 
modified  turbulence  model  in  low  Reynolds  number  flows,  in  which  the  original  model  overestimated  the 
length  of  transition.  The  location  of  transition  evidently  has  a  considerable  influence  on  computed  boundary 
layer  results.  Numerical  instabilities  had  to  be  overcome  partly  by  specifying  transition,  which  might  cause  a 
deviation  from  the  true  flow  behaviour  because  of  surpressed  laminar  separation  bubbles.  Since  the  method 
is  not  yet  capable  of  handling  highly  cambered  blade  sections  or  incidence  angles  corresponding  to  stall  and 
post-stall,  further  improvements  are  necessary  to  obtain  a  more  universal  code.  Future  work  should  attempt 
a  more  accurate  modelling  of  both  separating  and  transitional  flow  regimes.  In  the  near  future  progress  may 
be  made  by 

1.  Including  the  effects  of  viscous  wakes,  which  should  improve  our  ability  to  predict  separated  flow. 

2.  Relaxing  the  FLARE-approximation  by  introducing  an  upwind  differencing  scheme  for  streamwise 
derivatives. 

3.  Reconsidering  some  of  the  numerical  features,  like  interpolation  and  smoothing  algorithms. 

Ultimately,  a  satisfactory  analysis  of  viscous  cascade  flows  will  require 

4.  A  better  understanding  of  the  physics  of  transition,  eventually  leading  to  a  more  reliable  prediction  of 
transition. 

5.  A  more  advanced  interaction  model  to  resolve  the  difficulties  at  high  incidence  angles. 
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Fig.  5.  Lift  histories  of  a  converging  case.  Flow  past  a  staggered  cascade  at  high  Reynolds  numbe 
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Fig.  G.  Lift  histories  of  a  diverging  case.  Flow  past  a  staggered  cascade  at  low  Reynolds  number 
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Lift  histories  of  a  slightly  oscillating  case.   Flow  past  a  staggered  cascade 
at  low  Reynolds  number  with  a  changed  turbulence  model 
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Fig.    8.     Lift  histories  of  a  converging  case.    Flow  past  a  staggered  cascade  at  low 
Reynolds  number  with  a  changed  turbulence  model  and  fixed  transition 
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Fig.  10.  Development  of  pressure  distributions  with  number  of  iterations  for  a  converging  case 
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Fig.  11.  Development  of  local  skin  friction  distribution  with  number  of  iterations  for  a  converging  case 
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Fig.  12.  Development  of  local  skin  friction  distribution  with  number  of  iterations  for  a  com 


er£iii£  case 
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Fig.  13.  Development  of  displacement  thickness  distribution  with  number  of  iterations  for  a  diverging  case 
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Fij.  14.  Development  of  displacement  thickness  distribution  with  number  of  iterations  for  a  converging-  case 
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Fig.  15.  Lift  and  turning  angle  as  a  function  of  air  inlet  angle  for  two  different  airfoil  sections  in  cascade 
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Fig.  16.    Variation  of  pressure  distribution  with  air  inlet  angle  for  a  staggered  cascade 
at  high  Reynolds  number 


43 


"•> 

**£* 

^Zji 

^^.    Ijfl 

^v. 

""""•-•■ 

"--P"" 

^  J     t\]     I 

S*"v*-/jB»' 

'*■-.., 

** 

*•  i    v-  i^*^*! 

i^-ini 

"\.. 

N"Kv. 

^c^^ 

*r***- 

*■•*> 

"-•4  ^**i>^,-^ 

*., 

"""•■ 

-■■•IS  «S^     «J^ 

"••. 

^V  >^ 

<3>        ^ 

S    4 

<o 

& 

Fig.  17.     Variation  of  local  skin  friction  distribution  with  air  inlet  angle  for  a  stag- 
gered cascade  at  high  Reynolds  number 
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Fig.  18.    Variation  of  displacement  thickness  distribution  with  air  inlet  angle  for  a 
staggered  cascade  at  high  Reynolds  number 
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Cascade  Data: 

NACA  66-410 

Blada  Solidity:        0.78 

Stagger  Angia:  46  dag 


0.00  "■  ■  »  I  i  i  ■  I 

0.0    0.6       1.0  U/U 


Symbols: 


TR...  Transition 
SP...  Separation 


Reynolds  Number:  6000000 


Air  Inlet  Angle:     54.3  deg 


Fig.  19.    Velocity  profiles  oo  a  NACA  65-410  airfoil  in  cascade  at  high  Reynolds 
number  and  moderate  angle  of  attack 


Cascade  Data: 


Symbols: 


TR...  Transition 
SP...  Separation 


Reynolds  Number:  245000 


Air  Inlet  Angle:     54.3  deg 


Fig.  20.     Velocity  profiles  on  a  NACA  65-410  airfoil  in  cascade  at  low  Reynolds  num- 
ber and  moderate  anffle  of  attack 
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Scales 


Cascade  Data: 

NACA  66-410 

Blade  Solidity:        0.78 

Stagger  Anglo:  46  dag 


Symbols: 


TR...  Transition 
SP...  Separation 
RA...  Raattachmant 


Reynolds  Number  500000 


Air  Inlet  Angle:  48  deg 


Fig.  21.    Velocity  profiles  on  a  NACA  65-410  airfoil  in  cascade  at  low  Reynolds  num- 
ber and  low  angle  of  attack 


Cascade  Data: 

NACA  65-410 

Blade  Solidity:         0.78 

Staggar  Angle:  45  dag 


0.00 


0.0     0.5       1.0  U/U 


Symbols: 


TR...  Transition 
SP...  Separation 
RA...  Reattachment 


Air  Inlet  Angle:  58  deg 


Reynolds  Number:  500000 

Fig.  22.     Velocity  profiles  ou  a  NACA  65-410  airfoil  in  cascade  at  low  Reynolds  num- 
ber and  high  angle  of  attack 
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